diff --git a/src/database/data_container.hpp b/src/database/data_container.hpp index 6ecbff4..63b8b20 100644 --- a/src/database/data_container.hpp +++ b/src/database/data_container.hpp @@ -1,353 +1,355 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_DATABASE_DATACONTAINER_HPP #define SPECMICP_DATABASE_DATACONTAINER_HPP //! \file data_container.hpp //! \brief Storage class for the thermodynamics database #include "../types.hpp" #include "../physics/units.hpp" #include "species/component_list.hpp" #include "species/aqueous_list.hpp" #include "species/mineral_list.hpp" #include "species/gas_list.hpp" #include "species/sorbed_list.hpp" namespace specmicp { namespace database { extern const std::string water_label; //!< The water label in the database constexpr index_t water_id{0}; //!< Index of water in the database extern const std::string electron_label; //!< The electron label in the database constexpr index_t electron_id{1}; //!< Index of the electron in the database //! \brief Information used to trace the database struct Metadata { std::string path; //!< path to the database std::string name; //!< name of the database std::string version; //!< version of the database }; //! \brief Storage class - Contains the database //! //! The database is shared by (almost) everything else using a shared_ptr. //! This class is only a container. It should not be modified by hand but //! using the algorithms provided by specimicp::database::Database. //! //! \sa specimicp::database::Database struct DataContainer { - Metadata metadata; + DataContainer() {} + + Metadata metadata; //!< Metada, name path and version of the database //! \name Basis //! \brief The basis contains the components // ========================================== //! @{ ComponentList components; //!< The basis //! \brief Return the number of components in the basis index_t nb_component() const {return components.size();} //! \brief Return the number of all aqueous components, not including water and electron index_t nb_aqueous_components() const noexcept {return nb_component()-2;} //! \brief Return the id of component 'label' //! //! Return no_species if the component does not exist index_t get_id_component(const std::string& label) const { return components.get_id(label); } //! \brief Return the id of component 'label' std::string get_label_component(index_t id) const{ return components.get_label(id); } //! \brief Return the molar mass of 'component' in g/mol scalar_t molar_mass_basis(index_t component) const NOEXCEPT {return components.molar_mass(component);} //! \brief Return the molar mass of 'component' in kg/mol scalar_t molar_mass_basis_si(index_t component) const NOEXCEPT {return 1e-3*molar_mass_basis(component);} //! \brief Return the molar mass of 'component' in 'mass_unit'/mol scalar_t molar_mass_basis(index_t component, units::MassUnit mass_unit) const NOEXCEPT; //! \brief Return the charge of a component const scalar_t& charge_component(index_t i) const NOEXCEPT {return components.charge(i);} //! \brief Return the 'a_i' Debye-Huckel parameter for a component const scalar_t& a_debye_component(index_t i) const NOEXCEPT {return components.a_debye(i);} //! \brief Return the 'b_i' extended Debye-Huckel parameter for a component const scalar_t& b_debye_component(index_t i) const NOEXCEPT {return components.b_debye(i);} //! \brief Return the index of the water in the basis constexpr static index_t water_index() noexcept { return water_id; } //! \brief return the index of the electron in the basis constexpr static index_t electron_index() noexcept { return electron_id; } //! \brief Range over the components range_t range_component() const { return components.range();} //! \brief Range over the aqueous components range_t range_aqueous_component() const { return boost::irange(static_cast(2), nb_component());} //! @} //! \name Secondary aqueous species //! \brief The secondary aqueous species // ========================================= //! @{ AqueousList aqueous; //!< The list of aqueous species //! \brief Return the number of aqueous species in the system (not counting components) index_t nb_aqueous() const {return aqueous.size();} //! \brief Return the stoichiometric coefficient of 'i' in the dissociation reaction of species 'j' const scalar_t& nu_aqueous(index_t j, index_t i) const {return aqueous.nu_ji(j, i);} //! \brief Return the logk of dissociation reaction of species 'j' const scalar_t& logk_aqueous(index_t j) const {return aqueous.logk(j);} //! \brief Return the id of aqueous species 'label' //! //! Return no_species if the aqueous species does not exist index_t get_id_aqueous(const std::string& label) const{ return aqueous.get_id(label); } //! \brief Return the id of aqueous species 'label' std::string get_label_aqueous(index_t id) const { return aqueous.get_label(id); } //! \brief Return the charge of a secondary aqueous species const scalar_t& charge_aqueous(index_t j) const NOEXCEPT {return aqueous.charge(j);} //! \brief Return the 'a_i' Debye-Huckel parameter for a secondary aqueous species const scalar_t& a_debye_aqueous(index_t j) const NOEXCEPT {return aqueous.a_debye(j);} //! \brief Return the 'b_i' extended Debye-Huckel parameter for a secondary aqueous species const scalar_t& b_debye_aqueous(index_t j) const NOEXCEPT {return aqueous.b_debye(j);} //! \brief Return true if the corresponding equation is a half-cell reaction bool is_half_cell_reaction(index_t aqueous_species) const NOEXCEPT { return (nu_aqueous(aqueous_species, electron_index()) != 0.0); } //! \brief Range over the secondary aqueous species range_t range_aqueous() const { return aqueous.range();} //! @} //! \name Sorbed species //! \brief The adsorbed species // ========================================= //! @{ SorbedList sorbed; //!< The list of sorbed species //! \brief Return the number of sorbed species index_t nb_sorbed() const {return sorbed.size();} //! \brief Return the stoichiometric coefficient of 'i' in the dissociation reaction of species 'j' const scalar_t& nu_sorbed(index_t j, index_t i) const {return sorbed.nu_ji(j, i);} //! \brief Return the logk of dissociation reaction of species 'j' const scalar_t& logk_sorbed(index_t j) const {return sorbed.logk(j);} //! \brief Return the number of sorption sites occupied by a sorbed species const scalar_t& nb_sorption_sites(index_t sorbed_species) const { return sorbed.nb_sorption_site_occupied(sorbed_species); } //! \brief Return the id of sorbed species 'label' //! //! Return no_species if the sorbed species does not exist index_t get_id_sorbed(const std::string& label) const { return sorbed.get_id(label); } //! \brief Return the id of sorbed species 'label' std::string get_label_sorbed(index_t id) const { return sorbed.get_label(id); } //! \brief Range over the sorbed species range_t range_sorbed() const { return boost::irange((index_t) 0, nb_sorbed());} //! \brief Return true if the corresponding equation is a half-cell reaction bool is_sorbed_half_cell_reaction(index_t sorbed_species) const { return (nu_sorbed(sorbed_species, electron_index()) != 0.0); } //! @} //! \name Minerals //! \brief The solid phases // ========================================= //! @{ MineralList minerals; //!< The list of solid phases at equilibrium MineralList minerals_kinetic; //!< The list of solid phases governed by kinetic laws //! \brief Returns the number of solid phases at equilibrium index_t nb_mineral() {return minerals.size();} //! \brief Returns the number of solid phases governed by kinetic laws index_t nb_mineral_kinetic() {return minerals_kinetic.size();} //! \brief Return the stoichiometric coefficient for 'component' in the dissolution reaction of 'mineral' const scalar_t& nu_mineral(index_t mineral, index_t component) const { return minerals.nu_ji(mineral, component); } //! \brief Return the stoichiometric coefficient for 'component' in the dissolution reaction of 'mineral' const scalar_t& nu_mineral_kinetic(index_t mineral, index_t component) const { return minerals_kinetic.nu_ji(mineral, component); } //! \brief Return the logk for the dissolution reaction of 'mineral' const scalar_t& logk_mineral(index_t mineral) const { return minerals.logk(mineral); } //! \brief Return the logk for dissolution reaction of 'mineral' const scalar_t& logk_mineral_kinetic(index_t mineral) const { return minerals_kinetic.logk(mineral); } //! \brief Return true if the corresponding equation is a half-cell reaction bool is_mineral_half_cell_reaction(index_t mineral_species) const { return (nu_mineral(mineral_species, electron_index()) != 0.0); } //! \brief Return true if the corresponding equation is a half-cell reaction bool is_mineral_kinetic_half_cell_reaction(index_t mineral_species) const { return (nu_mineral_kinetic(mineral_species, electron_index()) != 0.0); } //! \brief Return the id of solid phase 'label' //! //! Return no_species if the solid phase does not exist index_t get_id_mineral(const std::string& label) const { return minerals.get_id(label); } //! \brief Return the id of solid phase 'label' std::string get_label_mineral(index_t id) const { return minerals.get_label(id); } //! \brief Return the id of solid phase 'label' //! //! Return no_species if the solid phase does not exist index_t get_id_mineral_kinetic(const std::string& label) const { return minerals_kinetic.get_id(label); } //! \brief Return the id of solid phase 'label' std::string get_label_mineral_kinetic(index_t id) const { return minerals_kinetic.get_label(id); } //! \brief Range over the solid phases (at equilibrium) range_t range_mineral() const { return minerals.range();} //! \brief Range over the solid phases governed by equilibrium range_t range_mineral_kinetic() const { return minerals_kinetic.range();} //! \brief Return the molar mass (kg/mol) of a mineral scalar_t molar_mass_mineral(index_t mineral) const; //! \brief Return the molar mass (kg/mol) of a mineral governed by kinetic scalar_t molar_mass_mineral_kinetic(index_t mineral) const; //! \brief Return the molar mass of 'mineral' in 'mass_uinit'/mol scalar_t molar_mass_mineral(index_t mineral, units::MassUnit mass_unit) const; //! \brief Return the molar mass of 'mineral_kinetic' in 'mass_uinit'/mol scalar_t molar_mass_mineral_kinetic(index_t mineral_kinetic, units::MassUnit mass_unit) const; //! \brief Return the scaling factor for the molar volume scalar_t scaling_molar_volume(units::LengthUnit length_unit) const; //! \brief Return the molar volume whatever it is scalar_t unsafe_molar_volume_mineral(index_t m) const { return minerals.molar_volume(m); } //! \brief Return the molar volume (m^3/mol) of a mineral scalar_t molar_volume_mineral(index_t m) const; //! \brief Return the molar volume (m^3/mol) of a mineral governed by kinetic scalar_t molar_volume_mineral_kinetic(index_t m) const; //! \brief Return the molar volume of a mineral in mol/(length_unit)^3 scalar_t molar_volume_mineral(index_t m, units::LengthUnit length_unit) const; //! @} //! \name Gas //! \brief The gas present in the system // ========================================= //! @{ GasList gas; //!< The list of gas //! \brief Return the number of gas in the system index_t nb_gas() const {return gas.size();} //! \brief Return the stoichiometric coefficient of 'component' in the dissolution reaction of 'gas_id' const scalar_t& nu_gas(index_t gas_id, index_t component) const {return gas.nu_ji(gas_id, component);} //! \brief Return the logk of he dissolution reaction of 'gas_id' const scalar_t& logk_gas(index_t gas_id) const {return gas.logk(gas_id);} //! \brief Return true if the corresponding equation is a half-cell reaction bool is_gas_half_cell_reaction(index_t gas_species) const { return (nu_gas(gas_species, electron_index()) != 0.0); } //! \brief Return the id of gas 'label' //! //! Return no_species if the gas does not exist index_t get_id_gas(const std::string& label) const { return gas.get_id(label); } //! \brief Return the id of solid phase 'label' std::string get_label_gas(index_t id) const { return gas.get_label(id); } //! \brief Range over the gas phases range_t range_gas() const { return boost::irange((index_t) 0, nb_gas());} //! @} //! \brief Return true if the database is equal to 'other' //! This is just a test on the names in the database not all the values bool operator ==(const DataContainer& other) { return (get_hash() == other.get_hash()); } // Status // ------ //! \brief Return true if the database is valid; bool is_valid() const; //! \brief Freeze the database void freeze_db() {m_fixed_hash = get_hash();} //! \brief Return a hash of the database //! This only take into account the labels, not the values size_t get_hash() const; private: size_t m_fixed_hash {0}; }; } // end namespace database } // end namespace specmicp #include namespace specmicp { namespace database { //! \brief Pointer to a database //! //! This is a smart pointer so we don't have to worry about memory leaks. //! It is intented to be shared between all the modules that needs it. using RawDatabasePtr = std::shared_ptr; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_DATACONTAINER_HPP diff --git a/src/database/reader.cpp b/src/database/reader.cpp index d893b34..1269171 100644 --- a/src/database/reader.cpp +++ b/src/database/reader.cpp @@ -1,650 +1,645 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 "reader.hpp" #include "json/json.h" #include #include #include #include #include #include #include "utils/log.hpp" #define INDB_SECTION_BASIS "Basis" #define INDB_SECTION_AQUEOUS "Aqueous" #define INDB_SECTION_MINERALS "Minerals" #define INDB_SECTION_GAS "Gas" #define INDB_SECTION_SORBED "Sorbed" #define INDB_ATTRIBUTE_LABEL "label" #define INDB_ATTRIBUTE_ACTIVITY "activity" #define INDB_ATTRIBUTE_ACTIVITY_A "a" #define INDB_ATTRIBUTE_ACTIVITY_B "b" #define INDB_ATTRIBUTE_MOLARMASS "molar_mass" #define INDB_ATTRIBUTE_MOLARVOLUME "molar_volume" #define INDB_ATTRIBUTE_COMPOSITION "composition" #define INDB_ATTRIBUTE_LOGK "log_k" #define INDB_ATTRIBUTE_NBSITEOCCUPIED "nb_sites_occupied" #define INDB_ATTRIBUTE_FLAG_KINETIC "flag_kinetic" #define INDB_OPEN_DELIMITER_CHARGE '[' #define INDB_CLOSE_DELIMITER_CHARGE ']' #define INDB_SECTION_METADATA "Metadata" #define INDB_ATTRIBUTE_NAME "name" #define INDB_ATTRIBUTE_VERSION "version" #define INDB_ELECTRON "E[-]" #include namespace specmicp { namespace database { - -static std::hash hash_fn; - //! \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); //! \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); // safe access to values std::string obtain_label(Json::Value& species); void obtain_composition(Json::Value& species, std::map& compo); scalar_t obtain_logk(Json::Value& species); bool is_kinetic(Json::Value& species); -// void DataReader::parse(const std::string& filepath) { std::ifstream datafile(filepath, std::ifstream::binary); if (not datafile.is_open()) { ERROR << "Database not found : " << filepath; std::string message("Database not found : " + filepath); throw std::invalid_argument(message); } data->metadata.path = filepath; - return parse(datafile); + parse(datafile); datafile.close(); } Json::Value DataReader::parse_json(std::istream& input) { Json::Reader reader; Json::Value root; DEBUG << "Parsing database"; bool parsingSuccessful = reader.parse(input, root); if ( !parsingSuccessful ) { // report to the user the failure and their locations in the document. std::string message("Failed to parse configuration\n" + reader.getFormattedErrorMessages()); throw std::runtime_error(message); } return root; } void DataReader::parse(std::istream& input ) { Json::Value root = parse_json(input); - m_electron_hash = hash_fn(INDB_ELECTRON); check_for_mandatory_value(root, INDB_SECTION_METADATA); parse_metadata(root[INDB_SECTION_METADATA]); // Basis // ------ check_for_mandatory_value(root, INDB_SECTION_BASIS); parse_basis(root[INDB_SECTION_BASIS]); // Aqueous species // --------------- check_for_mandatory_value(root, INDB_SECTION_AQUEOUS); AqueousList alist; parse_aqueous(root[INDB_SECTION_AQUEOUS], alist); data->aqueous = std::move(alist); // Minerals // -------- check_for_mandatory_value(root, INDB_SECTION_MINERALS); MineralList minerals, minerals_kinetic; parse_minerals(root[INDB_SECTION_MINERALS], minerals, minerals_kinetic); data->minerals = std::move(minerals); data->minerals_kinetic = std::move(minerals_kinetic); // Gas // ---- GasList glist(0, data->nb_component()); if (root.isMember(INDB_SECTION_GAS)) parse_gas(root[INDB_SECTION_GAS], glist); else glist.set_valid(); data->gas = std::move(glist); // Sorbed species // -------------- SorbedList slist(0, data->nb_component()); if (root.isMember(INDB_SECTION_SORBED)) parse_sorbed(root[INDB_SECTION_SORBED], slist); else slist.set_valid(); data->sorbed = std::move(slist); } void DataReader::parse_metadata(Json::Value &root) { check_for_mandatory_value(root, INDB_ATTRIBUTE_NAME); data->metadata.name = root[INDB_ATTRIBUTE_NAME].asString(); check_for_mandatory_value(root, INDB_ATTRIBUTE_VERSION); data->metadata.version = root[INDB_ATTRIBUTE_VERSION].asString(); } void DataReader::parse_basis(Json::Value& basis) { DEBUG << "Size basis : " << basis.size(); data->components = ComponentList(basis.size()); for (int id: data->components.range()) { Json::Value species = basis[id]; check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL); std::string label = obtain_label(species); if (id == water_id and label != water_label) { throw invalid_database("The first component should be "+water_label+"."); } else if (id == electron_id and label != electron_label) { throw invalid_database("The second component should be "+electron_label+"."); } auto is_already_in_the_database = data->components.get_id(label); if (is_already_in_the_database != no_species) throw invalid_database("Component : " + label + "is already in the database."); SPAM << "Parsing component : " << label; ComponentValues values; values.label = label; values.ionic_values.charge = charge_from_label(label); if (species.isMember(INDB_ATTRIBUTE_ACTIVITY)) { values.ionic_values.a_debye = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_A].asDouble(); values.ionic_values.b_debye = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_B].asDouble(); } check_for_mandatory_value(species, INDB_ATTRIBUTE_MOLARMASS); values.molar_mass = species[INDB_ATTRIBUTE_MOLARMASS].asDouble(); data->components.set_values(id, std::move(values)); } data->components.set_valid(); } void DataReader::parse_aqueous(Json::Value& aqueous, AqueousList& alist) { DEBUG << "Parse aqueous species"; alist = AqueousList(aqueous.size(), data->nb_component()); //const auto size = data->nb_aqueous(); for (auto id: alist.range()) { Json::Value& species = aqueous[static_cast(id)]; AqueousValues values; values.label = obtain_label(species); values.ionic_values.charge = charge_from_label(values.label); if (species.isMember(INDB_ATTRIBUTE_ACTIVITY)) { values.ionic_values.a_debye = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_A].asDouble(); values.ionic_values.b_debye = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_B].asDouble(); } values.logk = obtain_logk(species); alist.set_values(id, std::move(values)); // equation std::map compo; obtain_composition(species, compo); scalar_t test_charge {0.0}; std::map to_canonicalize; for (const auto& it: compo) { // component ? const auto search = data->components.get_id(it.first); if(search != no_species) { alist.set_nu_ji(id, search, it.second); test_charge += it.second*data->charge_component(search); } // aqueous species ? else { // in the document we parse ? const auto id_aq = alist.get_id(it.first); if (id_aq != no_species) { to_canonicalize.insert({id_aq, it.second}); test_charge += it.second*alist.charge(id_aq); } // cannot be found... else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + alist.get_label(id) +"."); } } } if (test_charge != alist.charge(id)) { throw db_invalid_data("Total charge is not zero for aqueous species : '" + alist.get_label(id) + "' - Charge-decomposition - charge species : "+ std::to_string(test_charge - alist.charge(id) )+"."); } for (const auto& tocanon: to_canonicalize) { alist.canonicalize(id, tocanon.first, tocanon.second); } } alist.set_valid(); } void DataReader::parse_minerals(Json::Value& minerals, MineralList& minerals_list, MineralList& minerals_kinetic_list) { DEBUG << "Parse minerals"; index_t nb_mineral = minerals.size(); // find kinetic minerals int nb_kin = 0; for (int id=0; idnb_component()); minerals_kinetic_list = MineralList(nb_kin, data->nb_component()); // id for each class of minerals index_t id_in_eq=0; index_t id_in_kin=0; for (index_t id=0; id(id)]; //check if material is at equilibrium or governed by kinetics bool is_kin = is_kinetic(species); MineralValue value; value.label = std::move(obtain_label(species)); value.logk = obtain_logk(species); if (species.isMember(INDB_ATTRIBUTE_MOLARVOLUME)) { value.molar_volume = species[INDB_ATTRIBUTE_MOLARVOLUME].asDouble(); } else { // ###FIXME value.molar_volume = -1; } auto& mlist = is_kin?minerals_kinetic_list:minerals_list; auto& true_id = is_kin?id_in_kin:id_in_eq; mlist.set_values(true_id, std::move(value)); // equation // -------- std::map compo; obtain_composition(species, compo); double test_charge = 0; std::map to_canonicalize; for (const auto& it: compo) { auto idc = data->components.get_id(it.first); if(idc != no_species) { // this is a component mlist.set_nu_ji(true_id, idc, it.second); test_charge += it.second * data->charge_component(idc); } else { // this is an aqueous species auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second * data->charge_aqueous(idaq); } else // or something we don't know { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + mlist.get_label(true_id) +"."); } } } if (test_charge != 0) { throw db_invalid_data("Total charge is not zero for mineral : '"+ mlist.get_label(true_id) + "' - total charge : "+std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { mlist.canonicalize(true_id, data->aqueous, tocanon.first, tocanon.second); } // update index ++true_id; } specmicp_assert(id_in_kin == minerals_kinetic_list.size()); specmicp_assert(id_in_eq == minerals_list.size()); minerals_list.set_valid(); minerals_kinetic_list.set_valid(); } void DataReader::parse_gas(Json::Value& gas, GasList& glist) { DEBUG << "Parse gas"; glist = GasList(gas.size(), data->nb_component()); for (auto id: glist.range()) { Json::Value& species = gas[static_cast(id)]; GasValues values; values.label = obtain_label(species); auto is_already_in_the_database = data->gas.get_id(values.label); if (is_already_in_the_database != no_species) throw invalid_database("Component : " + values.label + "is already in the database."); values.logk = obtain_logk(species); glist.set_values(id, std::move(values)); // equation std::map compo; obtain_composition(species, compo); scalar_t test_charge = 0.0; std::map to_canonicalize; for (auto it: compo) { const auto idc = data->components.get_id(it.first); // component if(idc != no_species) { glist.set_nu_ji(id, idc, it.second); test_charge += it.second*data->charge_component(idc); } // aqueous species else { const auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second*data->charge_aqueous(idaq); } else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + data->gas.get_label(id) +"."); } } } if (test_charge != 0) { throw db_invalid_data("Total charge is not zero for gas '"+data->gas.get_label(id)+ "' : " + std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { glist.canonicalize(id, data->aqueous, tocanon.first, tocanon.second); } } glist.set_valid(); } void DataReader::parse_sorbed(Json::Value& sorbed, SorbedList& slist) { DEBUG << "Parse sorbed species"; slist = SorbedList(sorbed.size(), data->nb_component()); for (auto id: slist.range()) { Json::Value& species = sorbed[static_cast(id)]; SorbedValues values; values.label = obtain_label(species); auto is_already_in_the_database = slist.get_id(values.label); if (is_already_in_the_database != no_species) throw invalid_database("Component : " + values.label + "is already in the database."); values.logk = obtain_logk(species); check_for_mandatory_value(species, INDB_ATTRIBUTE_NBSITEOCCUPIED); const auto nb_site_occupied = species[INDB_ATTRIBUTE_NBSITEOCCUPIED].asDouble(); if (nb_site_occupied < 0) { throw db_invalid_data("The number of sites occupied by a sorbed species must be positive. (Species : " + values.label + ", number of sites : " +std::to_string(nb_site_occupied) + ")"); } values.sorption_site_occupied = nb_site_occupied; slist.set_values(id, std::move(values)); std::map compo; obtain_composition(species, compo); double test_charge = 0; std::map to_canonicalize; for (auto it: compo) { const auto idc = data->components.get_id(it.first); if(idc != no_species) { slist.set_nu_ji(id, idc, it.second); test_charge += it.second*data->charge_component(idc); } else { const auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second*data->charge_aqueous(idaq); } else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + slist.get_label(id) +"."); } } } if (test_charge != 0) { throw db_invalid_data("Total charge is not zero for gas '"+slist.get_label(id)+ "' : " + std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { slist.canonicalize(id, data->aqueous, tocanon.first, tocanon.second); } } slist.set_valid(); } void parse_equation(const std::string& equation, std::map &compo) { std::vector list_compo; boost::split(list_compo, equation, [](char input){return input == ',';}, boost::token_compress_on); for (auto it=list_compo.begin(); it!=list_compo.end(); ++it) { std::string& toparse = *it; double coeff = 0; boost::trim_all(toparse); unsigned int pos_end = 0; while (pos_end < toparse.size()) { if (std::isalpha(toparse[pos_end])) // or std::isblank(toparse[pos_end])) { break; } ++pos_end; } std::string label(toparse.substr(pos_end, toparse.size()-pos_end)); boost::trim_all(label); if (pos_end == 0) coeff =1; else if (pos_end == 1 and toparse[0] == '-') coeff=-1; else { std::string tofloat = toparse.substr(0, pos_end); while (pos_end > 1) { if ((tofloat[0] == '-' or tofloat[0] == '+') and std::isspace(tofloat[1]) ) { tofloat = tofloat[0] + tofloat.substr(2,pos_end-2); --pos_end; continue; } else { break; } } if ( tofloat[pos_end-1] == '-' ) {coeff = -1;} else if (tofloat[pos_end-1] == '+') {coeff = +1;} else {coeff = std::stof(tofloat);} } compo[label] = coeff; } } double charge_from_label(const std::string& label) { int sign=1; auto start= label.end(); auto stop = label.end(); for (auto it=label.begin(); it != label.end(); ++it) { if (*it == INDB_OPEN_DELIMITER_CHARGE) { start = it; } } if (start == label.end()) {return 0;} // no charge specified for (auto it=start+1; it != label.end(); ++it) { if (*it == INDB_CLOSE_DELIMITER_CHARGE) { stop = it; } } if (stop == label.end()) {throw db_invalid_syntax("Charge : missing closing delimiter for species : "+label+".");} start = start+1; if (stop == start) {return 0;} // nothing inside bracket // get the signs if (*(stop-1) == '-') { sign = -1; stop =stop-1; } else if (*(stop-1) == '+') { sign = +1; stop = stop-1; } if (stop == start) { return sign; // just the sign, |z|=1 } double charge; try { charge = sign*std::stof(std::string(start,stop)); } catch (std::invalid_argument&) { throw db_invalid_syntax("Charge : bad formatting for species :"+ label+"."); } return charge; } 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 +"'."); } } inline Json::Value& get_mandatory_section(Json::Value& root, const std::string& section) { check_for_mandatory_value(root, section); return root[section]; } std::string obtain_label(Json::Value& species) { check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL); return species[INDB_ATTRIBUTE_LABEL].asString(); } void obtain_composition(Json::Value& species, std::map& compo) { check_for_mandatory_value(species, INDB_ATTRIBUTE_COMPOSITION); parse_equation(species[INDB_ATTRIBUTE_COMPOSITION].asString(), compo); } scalar_t obtain_logk(Json::Value& species) { check_for_mandatory_value(species, INDB_ATTRIBUTE_LOGK); return species[INDB_ATTRIBUTE_LOGK].asDouble(); } bool is_kinetic(Json::Value& species) { bool is_kin {false}; if (species.isMember(INDB_ATTRIBUTE_FLAG_KINETIC)) { is_kin = species[INDB_ATTRIBUTE_FLAG_KINETIC].asBool(); } return is_kin; } } // end namespace specmicp } // end namespace chemy diff --git a/src/database/reader.hpp b/src/database/reader.hpp index 9c39bf2..50210ed 100644 --- a/src/database/reader.hpp +++ b/src/database/reader.hpp @@ -1,141 +1,139 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 DATABASE_READER_H #define DATABASE_READER_H //! \file reader.hpp Read the database from a json file #include "json/json-forwards.h" #include "module.hpp" #include namespace specmicp { namespace database { //! \brief Read a json file containing a database class DataReader: public DatabaseModule { public: //! \brief Empty constructor DataReader() {} //! \brief Constructor DataReader(RawDatabasePtr data): DatabaseModule(data) {} //! \brief Constructor //! //! @param filepath string containing the path to the database DataReader(const std::string& filepath): DatabaseModule() { parse(filepath); } //! \brief Constructor //! //! @param input input stream that contains the database DataReader(std::istream& input): DatabaseModule() { parse(input); } //! Return the databes RawDatabasePtr get_database() {return data;} //! \brief Parse the basis section //! //! Contains the list of primary species void parse_basis(Json::Value& basis_root); //! \brief Parse the aqueous section //! //! Contains the list of secondary species void parse_aqueous(Json::Value& aqueous_root, AqueousList& alist); //! \brief Parse the mineral section //! //! Contains the list of minerals void parse_minerals( Json::Value& minerals, MineralList &minerals_list, MineralList &minerals_kinetic_list ); //! \brief Parse the gas section void parse_gas(Json::Value& gas_root, GasList &glist); //! \brief Parse the sorbed species section void parse_sorbed(Json::Value& sorbed_root, SorbedList &slist); //! \brief Parse a json input Json::Value parse_json(std::istream& input); private: //! \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(std::istream& input); //! \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(const std::string& filepath); //! \brief Parse the metadata section //! //! we don't do much with them for now.... void parse_metadata(Json::Value& root); - std::size_t m_electron_hash; - }; //! \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); } // end namespace database } // end namespace specmicp #endif // DATABASE_READER_H diff --git a/src/database/selector.hpp b/src/database/selector.hpp index bd3d729..9227ede 100644 --- a/src/database/selector.hpp +++ b/src/database/selector.hpp @@ -1,119 +1,119 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_DATABASE_SELECTOR_HPP #define SPECMICP_DATABASE_SELECTOR_HPP //! \file selector.hpp Select components in the database #include #include "module.hpp" namespace specmicp { namespace database { //! \brief Select components in the database //! //! This class reduce the database by removing some components //! and all the species associated with these components. class DatabaseSelector: public DatabaseModule { public: DatabaseSelector(RawDatabasePtr thedata): DatabaseModule(thedata), m_is_component_to_remove(thedata->nb_component(), false) {} //! \brief Remove compoment in the list "id_component_to_remove" //! //! \param id_component_to_remove list of id of the components to remove void remove_component(const std::vector& id_component_to_remove); //! \brief keep only components in the "id_to_keep" list //! //! \param id_to_keep list of id of the components to keep in the database void keep_only_component(const std::vector& id_to_keep); //! \brief Remove all gas phase void remove_all_gas(); //! \brief Remove all sorbed species void remove_all_sorbed(); //! \brief Remove some specific aqueous species //! //! This method should not be used in normal usage. //! Removing the aqueous species will break the consistency of the database //! //! \param id_aqueous list of id of the aqueous to remove void remove_aqueous(const std::vector& id_aqueous); private: //! \brief Analyse wich component should be removed void analyse_component(const std::vector& id_component_to_remove); void select_secondary( ReactiveSpeciesList* toselect, 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 the correct minerals governed by kinetics void select_minerals_kinetic(const std::vector& id_component_to_remove); //! \brief Select the gas phases void select_gas(const std::vector& id_component_to_remove); //! \brief Select the sorbed species void select_sorbed(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(index_t id) {return m_is_component_to_remove[id];} + bool is_component_to_remove(index_t id) const {return m_is_component_to_remove[id];} //! \brief Return the number of component to keep in the db - uindex_t nb_component_to_keep() {return m_nb_component_to_keep;} + uindex_t nb_component_to_keep() const {return m_nb_component_to_keep;} std::vector m_is_component_to_remove; // true if component will be removed uindex_t m_nb_component_to_keep; // total number of component after the purging is done }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_SELECTOR_HPP diff --git a/src/database/species/base_wrapper.cpp b/src/database/species/base_wrapper.cpp index f349832..21664ad 100644 --- a/src/database/species/base_wrapper.cpp +++ b/src/database/species/base_wrapper.cpp @@ -1,54 +1,54 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 "base_wrapper.hpp" namespace specmicp { namespace database { //! \brief replace imformation of 'new_ind' by those contained in 'old_ind' -void MatrixSpeciesWrapper::move_erase(index_t old_ind, index_t new_ind, const std::vector is_col_to_remove) { +void MatrixSpeciesWrapper::move_erase(index_t old_ind, index_t new_ind, const std::vector& is_col_to_remove) { specmicp_assert(static_cast(is_col_to_remove.size()) == m_matrix.cols()); index_t new_col {0}; for (index_t col=0; col, 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_DATABASE_BASEWRAPPER_HPP #define SPECMICP_DATABASE_BASEWRAPPER_HPP #include "../../types.hpp" //! \file base_wrapper.hpp //! //! \brief Base wrapper around vectors and matrices //! //! These classes encapsulate an Eigen Matrix (or vector) //! //! \defgroup database_species namespace specmicp { namespace database { //! \class MatrixSpeciesWrapper //! \brief A wrapper around a matrix, to be used by a SpeciesList (or a derived class) instance //! //! \ingroup database_species Describing the species in the database class MatrixSpeciesWrapper { public: //! \brief Default constructor MatrixSpeciesWrapper() {} //! \brief Constructor initializing the matrix MatrixSpeciesWrapper(index_t size, index_t nb_cols): m_matrix(Matrix::Zero(size, nb_cols)) {} //! \name Size // ============ //! \brief Return and set the size of the matrix //! @{ //! \brief Return the number of rows in the matrix index_t size() const {return m_matrix.rows();} //! \brief Return the number of columnns in the matrix index_t cols() const {return m_matrix.cols();} //! \brief Resize the matrix //! This is a resize only for the number of rows void resize(index_t size) { m_matrix.conservativeResize(size, Eigen::NoChange); } //! \brief Resize the matrix (both rows and colums) void resize(index_t rows, index_t cols) { m_matrix.conservativeResize(rows, cols); } // @} //! \name Getter // ============= //! \brief Return the values //! @{ //! \brief Return the element ['row', ['col'] const scalar_t& operator()(index_t row, index_t col) const { return m_matrix(row, col); } //! \brief Return a const reference to the matrix const Matrix& get_matrix() const { return m_matrix; } //! \brief Return the row of a matrix Eigen::MatrixBase::ConstRowXpr get_row(index_t k) const { return m_matrix.row(k); } //! \brief Return the element ['row','col'] const scalar_t& get_value(index_t row, index_t col) const { return m_matrix(row, col); } //! @} //! \name Setter // ============= //! \brief Set the values //! @{ //! \brief Add 'other' to row 'k' template void add_to_row(index_t k, const Eigen::MatrixBase& other) { m_matrix.row(k) += other; } //! \brief Multiply row 'k' by 'multiplier' template void multiply_by(const Eigen::MatrixBase& multiplier) { m_matrix = m_matrix*multiplier; } //! \brief Set the value of the matrix by copying 'init_matrix' template void set_matrix(const Eigen::MatrixBase& init_matrix) { m_matrix = init_matrix; } //! \brief Set the element ['row','col'] to 'value' void set_value(index_t row, index_t col, scalar_t value) { m_matrix(row, col) = value; } //! \brief Reset row 'k' void reset_row(index_t k) { m_matrix.row(k).setZero(); } //! @} //! \name Move // =========== //! \brief Move values inside, and outside, of the matrix //! @{ //! \brief replace imformation of 'new_ind' by those contained in 'old_ind' void move_erase(index_t old_ind, index_t new_ind) { m_matrix.row(new_ind) = m_matrix.row(old_ind); } //! \brief replace imformation of 'new_ind' by those contained in 'old_ind' - void move_erase(index_t old_ind, index_t new_ind, const std::vector is_col_to_remove); + void move_erase(index_t old_ind, index_t new_ind, const std::vector& is_col_to_remove); //! \brief Move the row to another matrix void move_erase_to(index_t ind, MatrixSpeciesWrapper& other, index_t other_ind) { other.m_matrix.row(other_ind) = m_matrix.row(ind); m_matrix.row(ind).setZero(); } //! \brief Swap two rows void swap(index_t ind, MatrixSpeciesWrapper& other, index_t other_ind) { m_matrix.row(ind).swap(other.m_matrix.row(other_ind)); } //! @} private: Matrix m_matrix; }; //! \class VectorSpeciesWrapper //! \brief A wrapper around a vector //! //! \ingroup database_species class VectorSpeciesWrapper { public: //! \brief Default constructor VectorSpeciesWrapper() {} //! \brief Initialise the vector //! \param size initial size of the vector VectorSpeciesWrapper(index_t size): m_vector(Vector::Zero(size)) {} //! \name Size // ============ //! \brief Return and set the size of the vector //! @{ //! \brief Return the size of the vector index_t size() const {return m_vector.rows();} //! \brief Resize the vector void resize(index_t size) { m_vector.conservativeResize(size); } //! @} //! \name Getter // ============= //! \brief Return the values //! @{ //! \brief return the value of index k const scalar_t& operator()(index_t k) const { return m_vector(k); } //! \brief Return the inner product of the vector and 'other' template scalar_t dot(const Eigen::MatrixBase& other) const { return m_vector.dot(other); } //! \brief Return the value of index k const scalar_t& get_value(index_t k) const { return m_vector(k); } //! \brief Return the vector const Vector& get_vector() const { return m_vector; } //! @} //! \name Setter // ============= //! \brief Set the values //! @{ //! \brief add 'vector' template void add(const Eigen::MatrixBase& vector) { m_vector += vector; } //! \brief Multiply the vector by 'multiplier' template void multiply_by(const Eigen::MatrixBase& multiplier) { m_vector = m_vector*multiplier; } //! \brief Set the value at index k void set_value(index_t k, scalar_t value) { m_vector(k) = value; } //! \brief Copy 'vector' to the vector template void set_vector(const Eigen::MatrixBase& vector) { m_vector = vector; } //! \brief Apply a linear transformation to the vector template void transform(const Eigen::MatrixBase& transform_matrix) { specmicp_assert(transform_matrix.rows() == transform_matrix.cols()); m_vector = transform_matrix*m_vector; } //! @} //! \name Move // =========== //! \brief Move values inside, and outside, of the matrix //! @{ //! \brief Move value of old_ind at new_ind void move_erase(index_t old_ind, index_t new_ind) { m_vector(new_ind) = m_vector(old_ind); } //! \brief Move row 'ind' to 'other' at row 'other_ind' void move_erase_to(index_t ind, VectorSpeciesWrapper& other, index_t other_ind) { other.m_vector(other_ind) = m_vector(ind); m_vector(ind) = 0; } //! @} private: Vector m_vector; }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_BASEWRAPPER_HPP diff --git a/src/dfpmsolver/driver.hpp b/src/dfpmsolver/driver.hpp index a1398a8..de52fa7 100644 --- a/src/dfpmsolver/driver.hpp +++ b/src/dfpmsolver/driver.hpp @@ -1,138 +1,138 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_DFPMSOLVER_DRIVER_HPP #define SPECMICP_DFPMSOLVER_DRIVER_HPP #include "../types.hpp" namespace specmicp { namespace dfpmsolver { //! \brief Base class for a driver //! //! Options should be a subclass of DriverOptions //! Performance should be a subclass of DriverPerformance //! template class Driver { public: Driver(Program& the_program): m_program(the_program), m_options(), m_scaling_is_initialized(false) { } Driver(Program& the_program, Options some_options): m_program(the_program), m_options(some_options), m_scaling_is_initialized(false) { } // basic program management // ------------------------ //! Return the program Program& program() {return m_program;} //! Return the number of equations - index_t get_neq() {return m_program.get_neq();} + index_t get_neq() const {return m_program.get_neq();} // common process // -------------- //! \brief rescale the update if needed //! //! Return the step length scalar_t is_step_too_long(Vector& update); // options // ------- //! \brief Return a read/write reference to the options Options& get_options() {return m_options;} //! \brief Return a read-only reference to the options const Options& get_options() const {return m_options;} // performance // ----------- //! \brief Return a const reference to the performance const Performance& get_perfs() const {return m_performance;} // Scaling // ------- void initialize_scaling() { if (not m_scaling_is_initialized) { m_scaling.resize(get_neq()); m_scaling.setOnes(); } } void set_scaling(const Vector& scale) { specmicp_assert(scale.rows() == get_neq()); m_scaling = scale; m_scaling_is_initialized = true; } const Vector& scaling() const {return m_scaling;} scalar_t scaling(index_t id_eq) const {return m_scaling(id_eq);} protected: //! \brief Read/write access to the performance Performance& get_perfs() {return m_performance;} private: Program& m_program; //!< The program Performance m_performance; //!< The performance Options m_options; //!< The options Vector m_scaling; //!< Scaling factor bool m_scaling_is_initialized; }; } // end namespace dfpmsolver } // end namespace specmicp // implementation // ============== #include "driver.inl" #endif // SPECMICP_DFPMSOLVER_DRIVER_HPP diff --git a/src/dfpmsolver/parabolic_driver.hpp b/src/dfpmsolver/parabolic_driver.hpp index b205aa0..ead108a 100644 --- a/src/dfpmsolver/parabolic_driver.hpp +++ b/src/dfpmsolver/parabolic_driver.hpp @@ -1,185 +1,185 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_DFPMSOLVER_PARABOLICDRIVER_HPP #define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP #include "../types.hpp" #include "../utils/sparse_solvers/sparse_solver.hpp" #include "driver.hpp" #include "parabolic_structs.hpp" namespace specmicp { namespace dfpmsolver { //! \brief The parabolic driver //! //! Parabolic driver for finite element programs template class ParabolicDriver: public Driver { public: using base=Driver; using base::program; using base::get_neq; using base::get_options; using base::get_perfs; using base::scaling; using base::initialize_scaling; ParabolicDriver(Program& the_program): Driver(the_program), m_velocity(Vector::Zero(the_program.get_tot_ndf())), m_solver(nullptr) {} //! \brief Solve a timestep of length dt ParabolicDriverReturnCode solve_timestep(scalar_t dt, Vector& displacement); //! \brief Restart the current timestep ParabolicDriverReturnCode restart_timestep(Vector& displacement); //! \brief Check if the solution has been found ParabolicDriverReturnCode check_convergence(); //! \brief Compute the residuals void compute_residuals(const Vector& displacement, Vector& residual) { compute_residuals(displacement, m_velocity, residual); } //! \brief Compute the residuals void compute_residuals(const Vector& displacement, const Vector& velocity, Vector& residual) { program().compute_residuals(displacement, velocity, residual); get_perfs().nb_call_residuals += 1; } //! \brief Compute the jacobian void compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian ); //! \brief Return the norm of the current residuals double norm_residuals() {return m_residuals.norm();} //! \brief Read/Write reference to the velocity vector const Vector& get_velocity() const {return m_velocity;} Vector& velocity() {return m_velocity;} void set_velocity(Vector& velocity_vector); //! \brief Reset the Sparse solver //! This function should be called when a new solver need to be used void reset_solver() {return m_solver.reset(nullptr);} //! \brief Call this function if the pattern of the jacobian has changed void jacobian_pattern_has_changed() {reset_solver();} //! \brief Initialize the computation void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement); //! \brief Strang Linesearch //! //! ref : //! - Matthies et al. (1979) //! - JHP course notes ParabolicLinesearchReturnCode strang_linesearch( Vector& update, Vector& displacements, scalar_t& lambda ); scalar_t residuals_0() {return m_norm_0;} private: void reset_velocity(); //! \brief Backtracking Linesearch //! //! ref : //! - Algo A6.3.1 : Dennis and Schnabel (1983) //! - Nocedal & Wrigth (2006) ParabolicLinesearchReturnCode backtracking_linesearch( Vector& update, Vector& displacements, scalar_t& lambda_out ); //! Update the variables in displacement void update_variable( const Vector& update, scalar_t lambda, Vector& displacement ); //! \brief Set the predictor void set_predictor(Vector& displacement); //! \brief perform the linesearch ParabolicDriverReturnCode linesearch( Vector &update, Vector &displacements ); //! Compute the residuals for the linesearch double compute_residuals_linesearch( Vector& update, scalar_t lambda, Vector& displacement ); double compute_residuals_strang_linesearch( Vector &update, scalar_t lambda, Vector &displacement ); scalar_t update_norm(const Vector& update); - double m_norm_0; - double m_current_dt; + double m_norm_0 {0.0}; + double m_current_dt {-1.0}; Eigen::VectorXd m_gradient; Eigen::VectorXd m_velocity; Eigen::VectorXd m_residuals; Eigen::VectorXd m_predictor; Eigen::SparseMatrix m_jacobian; sparse_solvers::SparseSolverPtr, Vector, Vector> m_solver; - bool m_velocity_is_initialized; + bool m_velocity_is_initialized {false}; }; } // end namespace dfpmsolver } // end namespace specmicp // implementation // ============== #include "parabolic_driver.inl" #endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP diff --git a/src/micpsolver/estimate_cond_number.hpp b/src/micpsolver/estimate_cond_number.hpp index 0cd19b4..5be5a38 100644 --- a/src/micpsolver/estimate_cond_number.hpp +++ b/src/micpsolver/estimate_cond_number.hpp @@ -1,117 +1,117 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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_MICPSOLVER_ESTIMATECONDNUMBER_HPP #define SPECMICP_MICPSOLVER_ESTIMATECONDNUMBER_HPP #include //! \file estimate_cond_number.hpp Estimate the condition number of a matrix namespace specmicp { //! \namespace micpsolver Namespace for the MiCP Solver namespace micpsolver { //! \brief Estimate the condition number of a dense square triangular matrix //! //! References : //! - \cite Cline1979 //! - \cite Hager1984 //! //! @param tmatrix a triangular view of a dense matrix //! @param maxit maximum number of iterations done by the algorithm (2 triangular solve by iteration) //! template double estimate_condition_number(Eigen::TriangularView tmatrix, int maxit=10) { const int n = tmatrix.cols(); Eigen::VectorXd x = 1/n*Eigen::VectorXd::Ones(n); Eigen::VectorXd y(n); int cnt = 0; y = tmatrix.solve(x); while (cnt < maxit) { for (int i=0; i= 0) y(i) = 1; else y(i) = -1; } tmatrix.solveInPlace(y); y.reverseInPlace(); // transpose int j; const double maxi = y.array().abs().maxCoeff(&j); const double gamma = y.dot(x); if (maxi <= gamma) break; // This is a local maximum x= Eigen::VectorXd::Zero(n); x(j) = 1; y = tmatrix.solve(x); ++cnt; } // norm tmatrix // ------------ double nnorm; if (mode == Eigen::Lower) { nnorm = std::abs(tmatrix(0, 0)); for (int i=1; i-1; --i) { double normrow = 0; for (int j=i; j(); } } // end namespace micpsolver } // end namespace specmicp #endif // SPECMICP_MICPSOLVER_ESTIMATECONDNUMBER_HPP diff --git a/src/micpsolver/micpsolver.hpp b/src/micpsolver/micpsolver.hpp index 8a8a113..e954733 100644 --- a/src/micpsolver/micpsolver.hpp +++ b/src/micpsolver/micpsolver.hpp @@ -1,241 +1,241 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 SPECMIC_MICPSOLVER_MICPSOLVER_HPP #define SPECMIC_MICPSOLVER_MICPSOLVER_HPP #include "../types.hpp" #include "micpsolver_base.hpp" #include "ncp_function.hpp" //! \file micpsolver.hpp //! \brief The MiCP solver namespace specmicp { namespace micpsolver { //! \class MiCPSolver //! \brief The MiCP Solver //! //! ### Mathematics //! //! Solve //! - \f$ G(u, v) = 0\f$ //! - \f$ 0 \leq v \perp H(u,v) \geq 0 \f$ //! //! Using the penalized Fisher-Burmeister C-function. //! This is a semismooth Newton solver with many features including : //! //! - Non-monotone linesearch //! - Scaling of the jacobian //! - Check of the condition number of the jacobian //! - Check of the descent directions //! - Crashing //! - ... //! //! References : //! - Munson et al (2001) \cite Munson2001 //! - Facchinei and Pang (2003) \cite Facchinei2003 //! //! ### Usage //! //! The options to customize the solver are contained in the struct //! specmicp::micpsolver::MiCPSolverOptions . //! //! \tparam program_t a subclass of MiCPProg //! //! To separate the program and the solver, the program is passed as a shared_ptr. //! The MiCPSolver is implemented using the CRTP pattern. //! The advantage is that we have compile time polymorphism instead of //! the runtime polymorphism when using the virtual functions. //! I am not sure of the effectiveness of this method (I have not run any benchmarks) //! but since the solver is destined to be called many times any small gain is good to have. //! Also, more compile-time optimization is possible (especially inlining). //! Anyway, that was fun to code... //! //! \sa specmicp::micpsolver::MiCPProgram template class MiCPSolver: public MiCPSolverBaseProgram { public: using base = MiCPSolverBaseProgram; using base::get_options; using base::get_program; using base::get_perfs; using base::get_neq; using base::get_neq_free; //! \brief Constructor //! //! \param prog smart pointer toward an instance of a program_t to solve MiCPSolver(std::shared_ptr prog): MiCPSolverBaseProgram(prog), m_residuals(prog->total_variables()), m_phi_residuals(Vector::Zero(prog->total_variables())), m_jacobian(prog->total_variables(),prog ->total_variables()) {} // Merit function // ############## // Residuals and jacobian // ---------------------- //! \brief Compute the residuals, use internal storage //! //! \param[in] x the variables void compute_residuals(const Vector& x) { base::compute_residuals(x, m_residuals); } //! \brief Compute the jacobian //! //! Assumes that the residual have been computed before //! //! \param[in] x the variables void compute_jacobian(Vector& x) { base::compute_jacobian(x, m_jacobian); } //! \brief Reformulation of the residuals //! //! Reformulate the problem - assumes that r, the residual, has been computed before //! //! \param[in] x the variables //! \param[in] r the residuals //! \param[out] r_phi a vector of size neq, which will contain the reformulated residuals void reformulate_residuals(const Vector& x, const Vector& r, Vector& r_phi); //! \brief Reformulation of the residuals *inplace* //! //! \param[in] x the variables //! \param[in,out] r the residual, will contain the reformulated residuals void reformulate_residuals_inplace(const Vector &x, Vector &r); //! \brief Reformulation of the jacobian //! //! r is the original vector of residuals (not reformulated) void reformulate_jacobian(const Vector& x) { base::reformulate_jacobian_cck(x, m_residuals, m_jacobian); } // Algorithm // ######### //! \brief Solver the program using x as initial guess //! //! \param[in,out] x the initial guess, as output contains the solution (from the last iteration) MiCPSolverReturnCode solve(Vector& x); //! \brief Setup the residuals //! //! \param[in] x the current solution void setup_residuals(const Eigen::VectorXd& x) { compute_residuals(x); reformulate_residuals(x, m_residuals, m_phi_residuals); } //! \brief Setup the jacobian //! //! \param[in] x the current solution void setup_jacobian(Eigen::VectorXd& x) { compute_jacobian(x); reformulate_jacobian(x); m_grad_phi = m_jacobian.transpose() * m_phi_residuals; } //! \brief Solve the Newton system //! //! The system will be scaled before being solved, may be expensive and not necessary. //! Do disable the scaling, disable the corresponding options in specmicp::MiCPSolver::MiCPSolverOptions. //! It assumes that the Newton system has been formed //! //! \param[out] update the update to apply to the solution //! //! \sa search_direction_calculation_no_scaling. MiCPSolverReturnCode search_direction_calculation(Vector& update); //! \brief Solve the Newton system - does not scale the jacobian //! //! The system will not be scaled. //! It assumes that the Newton system has been formed //! //! \param[out] update the update to apply to the solution //! //! \sa search_direction_scaling MiCPSolverReturnCode search_direction_calculation_no_scaling(Vector &update); //! \brief Linesearch //! //! It is a backtracking linesearch. //! If the correct option is set, this is a nonmonotone linesearch. //! //! \param[in,out] update the update of the timestep //! \param[in,out] x the current solution //! //! References : //! - \cite Dennis1983 //! - \cite Munson2001 int linesearch(Eigen::VectorXd& update, Eigen::VectorXd& x); //! \brief Crashing //! //! This function improves, if possible, the initial guess //! //! \param[in] x the current solution //! //! Reference : //! - \cite Munson2001 void crashing(Vector &x); private: // Residuals and jacobian Eigen::VectorXd m_residuals; //!< The residuals Eigen::VectorXd m_phi_residuals; //!< The reformulated residuals Eigen::VectorXd m_grad_phi; //!< The gradient of the reformulated residuals Eigen::MatrixXd m_jacobian; //!< The jacobian std::vector m_max_merits; //!< Contains the m best value of the merit function - bool m_gradient_step_taken; //!< True if the update was computed using the gradient + bool m_gradient_step_taken {false}; //!< True if the update was computed using the gradient }; } // end namespace micpsolver } // end namespace specmicp // ###############// // Implementation // // ###############// #include "micpsolver.inl" #endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP diff --git a/src/odeint/runge_kutta_step.inl b/src/odeint/runge_kutta_step.inl index 02b5a61..b3043a7 100644 --- a/src/odeint/runge_kutta_step.inl +++ b/src/odeint/runge_kutta_step.inl @@ -1,125 +1,124 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 #include "runge_kutta_step.hpp" // syntaxic coloration only #endif namespace specmicp { namespace odeint { template void EmbeddedRungeKuttaStep45::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, 3)*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); if (butcher_tableau_t::RK_evaluations == 6) { 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 ); } else { vec_t ak7; ytemp = y + timestep*(butcher().b(7,1)*dydx + butcher().b(7, 2)*ak2 + butcher().b(7, 3)*ak3 + butcher().b(7, 4)*ak4 + butcher().b(7, 5)*ak5 + butcher().b(7, 6)*ak6); get_rhs(x + butcher().a(7)*timestep, ytemp, ak7); 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 + butcher().c(7)*ak7); 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 + (butcher().c(7) - butcher().cs(7))*ak7 ); } if (get_options().non_negativity) set_non_negative(y_out); } template void EmbeddedRungeKuttaStep45::rk_step(vector_t& y, const vector_t& dydx, double& x, StepLength& stepl) { vector_t y_err; vector_t y_temp; double h = stepl.test; - double xnew; scalar_t errmax; for (int i=0; i(y_err) / get_options().eps; if (errmax < 1) break; double htemp = get_options().safety*h*std::pow(errmax, get_options().p_shrink); h = (h >= 0.0 ? std::max(htemp, 0.1*h) : std::min(htemp, 0.1*h)); - xnew = x+h; + double xnew = x+h; if (xnew == x) {throw std::runtime_error("stepsize underflow");} } stepl.did = h; x += h; if (errmax > get_options().errcon) {stepl.next = get_options().safety*h*std::pow(errmax, get_options().p_grow);} else {stepl.next = 5.0*h;} y = y_temp; } } // end namespace odeint } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index 5fe20ef..8c6bca3 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,1295 +1,1295 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. 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: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. 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. 3. Neither the name of the copyright holder 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 HOLDER 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 #include "adimensional_system.hpp" #include "../../utils/log.hpp" #include "../../physics/constants.hpp" #include "../../physics/laws.hpp" #include "adimensional_system_solution.hpp" #include #include // uncomment to activate the finite difference jacobian // #define SPECMICP_DEBUG_EQUATION_FD_JACOBIAN namespace specmicp { constexpr scalar_t log10 = std::log(10.0); // Constructor // =========== // No previous solution // -------------------- AdimensionalSystem::AdimensionalSystem( RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemOptions& options, const units::UnitsSet& units_set ): AdimemsionalSystemNumbering(ptrdata), OptionsHandler(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), m_second(ptrdata), m_equations(total_dofs(), ptrdata) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+1); number_eq(constraints); } // Previous solution // ----------------- AdimensionalSystem::AdimensionalSystem( RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemOptions& options, const units::UnitsSet& units_set ): AdimemsionalSystemNumbering(ptrdata), OptionsHandler(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), m_second(previous_solution), m_equations(total_dofs(), ptrdata) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+1); number_eq(constraints); } // Secondary variables constructor // =============================== // No previous solution // -------------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( const RawDatabasePtr& data ): secondary_molalities(Vector::Zero(data->nb_aqueous())), loggamma(Vector::Zero(data->nb_component()+data->nb_aqueous())), gas_fugacity(Vector::Zero(data->nb_gas())), gas_concentration(Vector::Zero(data->nb_gas())), sorbed_concentrations(Vector::Zero(data->nb_sorbed())) {} // Previous solution // ----------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( const AdimensionalSystemSolution& previous_solution ): secondary_molalities(previous_solution.secondary_molalities), loggamma(previous_solution.log_gamma), gas_fugacity(previous_solution.gas_fugacities), gas_concentration(Vector::Zero(previous_solution.gas_fugacities.rows())), sorbed_concentrations(previous_solution.sorbed_molalities) {} // IdEquations constructor // ======================= AdimensionalSystem::IdEquations::IdEquations( index_t nb_dofs, const RawDatabasePtr& data ): ideq(nb_dofs, no_equation), component_equation_type(data->nb_component()+1, no_equation), fixed_activity_species(data->nb_component()+1, no_species), active_aqueous(data->nb_aqueous(), false), active_gas(data->nb_gas(), false), active_sorbed(data->nb_sorbed()) {} // Equation numbering // ================== void AdimensionalSystem::number_eq( const AdimensionalSystemConstraints& constraints ) { index_t neq = 0; // Water // ===== if (constraints.water_equation != WaterEquationType::NoEquation) { m_equations.type_equation(dof_water()) = static_cast(constraints.water_equation); if (constraints.water_equation == WaterEquationType::MassConservation) { m_fixed_values(dof_water()) = constraints.total_concentrations(dof_water()); } m_equations.add_equation(dof_water(), &neq); } // Aqueous components // ================== number_eq_aqueous_component(constraints, neq); // Surface model // ============= if (constraints.surface_model.model_type == SurfaceEquationType::Equilibrium) { // add the equation m_equations.add_equation(dof_surface(), &neq); m_equations.type_equation(dof_surface()) = static_cast(constraints.surface_model.model_type); // setup the total concentration m_fixed_values(dof_surface()) = constraints.surface_model.concentration; } // Secondary species // ================= // Secondary aqueous species // ------------------------- bool include_half_cell_reaction = (constraints.electron_constraint.equation_type != ElectronEquationType::NoEquation); bool solve_electron_equation {false}; for (auto j: m_data->range_aqueous()) { bool can_exist { true }; if ( include_half_cell_reaction or not m_data->is_half_cell_reaction(j)) for (const auto& k: m_equations.nonactive_component) { if (m_data->nu_aqueous(j, k) != 0.0) { can_exist = false; break; } } else { can_exist = false; } m_equations.set_aqueous_active(j, can_exist); if (can_exist and m_data->is_half_cell_reaction(j)) solve_electron_equation = true; //std::cout << m_data->labels_aqueous[j] << "can exist ? " << can_exist << " - logk : " << m_data->logk_aqueous(j) << std::endl; } // Gas // --- for (index_t k: m_data->range_gas()) { bool can_exist = true; for (const index_t& n: m_equations.nonactive_component) { if (m_data->nu_gas(k, n) != 0.0) { can_exist = false; break; } } m_equations.set_gas_active(k, can_exist); } // Sorbed species // -------------- for (index_t s: m_data->range_sorbed()) { // Check if the surface model is computed if (constraints.surface_model.model_type != SurfaceEquationType::Equilibrium) { m_equations.set_sorbed_active(s, false); continue; } // If so, check that all components of the sorbed species exist bool can_exist = true; for (const index_t& k: m_equations.nonactive_component) { if (m_data->nu_sorbed(s, k) != 0.0) { can_exist = false; break; } } m_equations.set_sorbed_active(s, can_exist); } // Electron equation // ----------------- if (solve_electron_equation) { m_equations.add_equation(dof_electron(), &neq); m_equations.type_equation(dof_electron()) = static_cast(constraints.electron_constraint.equation_type); if (constraints.electron_constraint.equation_type == ElectronEquationType::Equilibrium) { m_fixed_values(dof_electron()) = 0.0; } else if (constraints.electron_constraint.equation_type == ElectronEquationType::FixedpE) { m_fixed_values(dof_electron()) = constraints.electron_constraint.fixed_value; m_equations.related_species(dof_electron()) = constraints.electron_constraint.species; //assert(m_fixed_activity_species[dof_electron()] >= 0 // and m_fixed_activity_species[dof_electron()] < m_data->nb_aqueous()); //assert(m_data->is_half_cell_reaction(m_fixed_activity_species[dof_electron()])); } // scaling if (get_options().scaling_electron == 0.0) { for (index_t component : m_data->range_aqueous_component()) { if (aqueous_component_equation_type(component) == AqueousComponentEquationType::MassConservation) { get_options().scaling_electron = total_concentration_bc(component); break; } } } } // above equations are 'free' (i.e. non constrained) m_equations.nb_free_variables = neq; // following equations are complementarity conditions // Minerals // ======== m_scaling_molar_volume = m_data->scaling_molar_volume(get_units().length); for (index_t m: m_data->range_mineral()) { bool can_precipitate = true; // just check that the molar volume exist auto molar_volume = m_data->molar_volume_mineral(m); // Remove minerals that cannot precipitate for (index_t& k: m_equations.nonactive_component) { if (m_data->nu_mineral(m, k) != 0.0 and molar_volume > 0.0) { can_precipitate = false; break; // this is not a mineral that can precipitate } } if (can_precipitate) { m_equations.add_equation(dof_mineral(m), &neq); } } m_equations.nb_tot_variables = neq; m_equations.nb_complementarity_variables = m_equations.nb_tot_variables - m_equations.nb_free_variables; } void AdimensionalSystem::number_eq_aqueous_component( const AdimensionalSystemConstraints& constraints, index_t& neq ) { using EqT = AqueousComponentEquationType; // First set the charge keeper if (constraints.charge_keeper != no_species) { if (constraints.charge_keeper == 0 or constraints.charge_keeper > m_data->nb_component()) { throw std::invalid_argument("The charge keeper must be an aqueous component. Invalid argument : " + std::to_string(constraints.charge_keeper)); } m_equations.type_equation(dof_component(constraints.charge_keeper)) = static_cast(EqT::ChargeBalance); } // Then go over fix fugacity gas for (const auto& it: constraints.fixed_fugacity_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component) + "' is already constrained, a fixed fugacity condition can not be applied"); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedFugacity); m_fixed_values(it.id_component) = it.log_value; m_equations.related_species(it.id_component) = it.id_gas; } // Then over the fix activity species for (const auto& it: constraints.fixed_activity_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component) + "' is already constrained, a fixed activity condition can not be applied."); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedActivity); m_fixed_values(it.id_component) = it.log_value; } // Finally number the equations for (index_t component: m_data->range_aqueous_component()) { // If no equation is assigned yet if (m_equations.type_equation(dof_component(component)) == static_cast(EqT::NoEquation)) { // Mass is conserved for this component //###FIXME: H[+], HO[-] const scalar_t& total_concentration = constraints.total_concentrations(dof_component(component)); if (std::abs(total_concentration) > get_options().cutoff_total_concentration) { m_equations.type_equation(dof_component(component)) = static_cast(EqT::MassConservation); m_fixed_values(dof_component(component)) = total_concentration; m_equations.add_equation(component, &neq); } else // add component to the nonactive component list { m_equations.add_non_active_component(component); } } // If equation is already assigned else { m_equations.add_equation(component, &neq); } } if (stdlog::ReportLevel() >= logger::Debug and m_equations.nonactive_component.size() > 0) { // if in debug mode list the non active components DEBUG << "Non active components :"; for (auto it: m_equations.nonactive_component) { DEBUG << " - " << it; } } } // ================ // // // // Residuals // // // // ================ // scalar_t AdimensionalSystem::weigthed_sum_aqueous(index_t component) const { scalar_t sum = 0.0; for (index_t aqueous: m_data->range_aqueous()) { if (not is_aqueous_active(aqueous)) continue; sum += m_data->nu_aqueous(aqueous,component)*secondary_molality(aqueous); } return sum; } scalar_t AdimensionalSystem::diff_weigthed_sum_aqueous(index_t diff_component, index_t component) const { scalar_t sum = 0.0; for (index_t aqueous: m_data->range_aqueous()) { if (not is_aqueous_active(aqueous)) continue; sum += log10*m_data->nu_aqueous(aqueous,diff_component)*m_data->nu_aqueous(aqueous,component)*secondary_molality(aqueous); } return sum; } scalar_t AdimensionalSystem::weigthed_sum_sorbed(index_t component) const { scalar_t sum = 0.0; for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; sum += m_data->nu_sorbed(s, component)*sorbed_species_concentration(s); } return sum; } scalar_t AdimensionalSystem::diff_weigthed_sum_sorbed(index_t diff_component, index_t component) const { scalar_t sum = 0.0; for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; sum += log10*m_data->nu_sorbed(s, diff_component)*m_data->nu_sorbed(s, component)*sorbed_species_concentration(s); } return sum; } scalar_t AdimensionalSystem::diff_surface_weigthed_sum_sorbed(index_t component) const { scalar_t sum = 0.0; for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; sum += log10*m_data->nb_sorption_sites(s)*m_data->nu_sorbed(s, component)*sorbed_species_concentration(s); } return sum; } scalar_t AdimensionalSystem::weigthed_sum_mineral(const Vector& x, index_t component) const { scalar_t sum = 0.0; for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation or m_data->nu_mineral(m, component) == 0.0) continue; const auto concentration = saturation_mineral(x, m)/molar_volume_mineral(m); sum += m_data->nu_mineral(m, component)*concentration; } return sum; } scalar_t AdimensionalSystem::weigthed_sum_gas(index_t component) const { scalar_t sum = 0.0; for (index_t k: m_data->range_gas()) { if (not is_active_gas(k) or m_data->nu_gas(k, component) == 0.0) continue; sum += m_data->nu_gas(k, component)*gas_concentration(k); } return sum; } scalar_t AdimensionalSystem::diff_weigthed_sum_gas(index_t diff_component, index_t component) const { scalar_t sum = 0.0; for (index_t k: m_data->range_gas()) { if (not is_active_gas(k) or m_data->nu_gas(k, component) == 0.0) continue; sum += log10*m_data->nu_gas(k, diff_component)*m_data->nu_gas(k, component)*gas_concentration(k); } return sum; } scalar_t AdimensionalSystem::residual_water(const Vector& x) const { scalar_t res {0.0}; if (water_equation_type() == WaterEquationType::MassConservation) { const scalar_t conc_w = density_water()*saturation_water(x); res = total_concentration_bc(0); res -= conc_w/m_data->molar_mass_basis_si(0); res -= conc_w*weigthed_sum_aqueous(0); if (ideq_surf() != no_equation) res -= conc_w*weigthed_sum_sorbed(0); res -= weigthed_sum_mineral(x, 0); if (m_data->nb_gas() > 0) res -= weigthed_sum_gas(0); res /= total_concentration_bc(0); } else if (water_equation_type() == WaterEquationType::SaturatedSystem) { res = 1 - saturation_water(x) - m_inert_volume_fraction; for (index_t mineral: m_data->range_mineral()) { res -= saturation_mineral(x, mineral); } } return res; } scalar_t AdimensionalSystem::residual_component(const Vector &x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::MassConservation); const scalar_t conc_w = density_water()*saturation_water(x); scalar_t res = total_concentration_bc(component); res -= conc_w*component_molality(x, component); res -= conc_w*weigthed_sum_aqueous(component); if (ideq_surf() != no_equation) res -= conc_w*weigthed_sum_sorbed(component); res -= weigthed_sum_mineral(x, component); if (m_data->nb_gas() > 0) res -= weigthed_sum_gas(component); res /= total_concentration_bc(component); return res; } scalar_t AdimensionalSystem::residual_fixed_activity(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedActivity); scalar_t res = (fixed_activity_bc(component) - log_gamma_component(component) - log_component_molality(x, component) ); res /= fixed_activity_bc(component); return res; } scalar_t AdimensionalSystem::residual_fixed_fugacity(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedFugacity); index_t id_g = m_equations.fixed_activity_species[component]; scalar_t res = fixed_fugacity_bc(component) + m_data->logk_gas(id_g); for (index_t component: m_data->range_aqueous_component()) { if (m_data->nu_gas(id_g, component) == 0) continue; res -= m_data->nu_gas(id_g, component)*( log_gamma_component(component) + log_component_molality(x, component)); } res /= fixed_fugacity_bc(component); return res; } scalar_t AdimensionalSystem::residual_mineral(const Vector& x, index_t m) const { specmicp_assert(ideq_min(m) != no_equation); scalar_t res = m_data->logk_mineral(m); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_mineral(m, i) != 0) { const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i); res -= m_data->nu_mineral(m, i)*log_activity_i; } } if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m)) res -= m_data->nu_mineral(m, m_data->electron_index())*log_activity_electron(x); return res; } scalar_t AdimensionalSystem::residual_charge_conservation(const Vector& x) const { scalar_t res = 0.0; for (index_t i: m_data->range_aqueous_component()) { if (m_data->charge_component(i) != 0 and ideq_paq(i) != no_equation) res += m_data->charge_component(i)*component_molality(x, i); } for (index_t j: m_data->range_aqueous()) { if (m_data->charge_aqueous(j) == 0 and not is_aqueous_active(j)) continue; res += m_data->charge_aqueous(j)*secondary_molality(j); } return res; } scalar_t AdimensionalSystem::residual_surface(const Vector &x) const { specmicp_assert(ideq_surf() != no_equation); const scalar_t conc_w = density_water()*saturation_water(x); scalar_t res = surface_total_concentration() - conc_w*free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; res -= conc_w*m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } return res/surface_total_concentration(); } scalar_t AdimensionalSystem::residual_electron(const Vector &x) const { specmicp_assert(electron_equation_type() == ElectronEquationType::Equilibrium); const scalar_t conc_w = density_water()*saturation_water(x); scalar_t res = 0.0; res -= conc_w * weigthed_sum_aqueous(m_data->electron_index()); if (ideq_surf() != no_equation) res -= conc_w * weigthed_sum_sorbed(m_data->electron_index()); res -= weigthed_sum_mineral(x, m_data->electron_index()); if (m_data->nb_gas() > 0) res -= weigthed_sum_gas(m_data->electron_index()); return res/get_options().scaling_electron; } void AdimensionalSystem::get_residuals(const Vector& x, Vector& residual) { residual.resize(total_variables()); // initialisation of 'main' secondary variables // They are especially nessary if the finite difference jacobian is used // and for the linesearch set_secondary_concentration(x); set_sorbed_concentrations(x); // // water if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x); // aqueous component for (index_t i: m_data->range_aqueous_component()) { switch (aqueous_component_equation_type(i)) { case AqueousComponentEquationType::NoEquation: break; case AqueousComponentEquationType::MassConservation: residual(ideq_paq(i)) = residual_component(x, i); break; case AqueousComponentEquationType::ChargeBalance: residual(ideq_paq(i)) = residual_charge_conservation(x); break; case AqueousComponentEquationType::FixedActivity: residual(ideq_paq(i)) = residual_fixed_activity(x, i); break; case AqueousComponentEquationType::FixedFugacity: residual(ideq_paq(i)) = residual_fixed_fugacity(x, i); break; } } // surface if (ideq_surf() != no_equation) residual(ideq_surf()) = residual_surface(x); // mineral for (index_t m: m_data->range_mineral()) { if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m); } // electron if (ideq_electron() != no_equation ) residual(ideq_electron()) = residual_electron(x); // std::cout << residual << std::endl; } // ================ // // // // Jacobian // // // // ================ // void AdimensionalSystem::get_jacobian(Vector& x, Matrix& jacobian) // //non-optimized Finite difference, for test only #ifdef SPECMICP_DEBUG_EQUATION_FD_JACOBIAN { finite_difference_jacobian(x, jacobian); return; } #else // analytical jacobian { analytical_jacobian(x, jacobian); return; } #endif void AdimensionalSystem::finite_difference_jacobian(Vector& x, Matrix& jacobian) { const int neq = total_variables(); Eigen::VectorXd res(neq); Eigen::VectorXd perturbed_res(neq); jacobian.setZero(neq, neq); get_residuals(x, res); for (int j=0; jmolar_mass_basis_si(0); tmp -= weigthed_sum_aqueous(0); tmp -= weigthed_sum_sorbed(0); tmp *= rho_w; jacobian(idw, idw) = tmp/factor; const scalar_t conc_w = density_water()*saturation_water(x); for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp = 0.0; tmp -= diff_weigthed_sum_aqueous(k, 0); tmp -= diff_weigthed_sum_sorbed(k, 0); // fixme gas tmp *= conc_w; tmp -= diff_weigthed_sum_gas(k, 0); jacobian(idw, ideq_paq(k)) = tmp/factor; } if (ideq_electron() != no_equation) { scalar_t tmp = 0.0; tmp -= diff_weigthed_sum_aqueous(m_data->electron_index(), 0); tmp -= diff_weigthed_sum_sorbed(m_data->electron_index(), 0); tmp*= conc_w; tmp -= diff_weigthed_sum_gas(m_data->electron_index(), 0); jacobian(idw, ideq_electron()) = tmp/factor; } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idw, ideq_min(m)) = -m_data->nu_mineral(m, 0)/molar_volume_mineral(m)/factor; } } else if (water_equation_type() == WaterEquationType::SaturatedSystem) { const index_t idw = ideq_w(); jacobian(idw, idw) = -1; for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idw, ideq_min(m)) = -1; } } } void AdimensionalSystem::jacobian_aqueous_components(Vector& x, Matrix& jacobian) { for (index_t i: m_data->range_aqueous_component()) { const index_t idp = ideq_paq(i); if (idp == no_equation) continue; switch (aqueous_component_equation_type(i)) { case AqueousComponentEquationType::NoEquation: continue; // Mass balance equation // ===================== case AqueousComponentEquationType::MassConservation: { const scalar_t conc_w = density_water()*saturation_water(x); const scalar_t factor = total_concentration_bc(i); // Aqueous components for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp_iip = 0; if (k == i) tmp_iip -= component_molality(x, i)*log10; tmp_iip -= diff_weigthed_sum_aqueous(k, i); tmp_iip -= diff_weigthed_sum_sorbed(k, i); tmp_iip *= conc_w; tmp_iip -= diff_weigthed_sum_gas(k, i); jacobian(idp, ideq_paq(k)) = tmp_iip/factor; } // Minerals for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idp, ideq_min(m)) = - m_data->nu_mineral(m, i)/molar_volume_mineral(m)/factor; } // Water if (ideq_w() != no_equation) { scalar_t tmp_iw = -component_molality(x, i); tmp_iw -= weigthed_sum_aqueous(i); tmp_iw -= weigthed_sum_sorbed(i); tmp_iw *= density_water(); jacobian(idp, ideq_w()) = tmp_iw/factor; } // Surface if (ideq_surf() != no_equation) { scalar_t tmp_s = -conc_w*diff_surface_weigthed_sum_sorbed(i); jacobian(idp, ideq_surf()) = tmp_s/factor; } // Electron if (ideq_electron() != no_equation) { scalar_t tmp = 0.0; tmp -= diff_weigthed_sum_aqueous(m_data->electron_index(), i); tmp -= diff_weigthed_sum_sorbed(m_data->electron_index(), i); tmp*= conc_w; tmp -= diff_weigthed_sum_gas(m_data->electron_index(), i); jacobian(idp, ideq_electron()) = tmp/factor; } break; } // Charge balance equation // ======================= case AqueousComponentEquationType::ChargeBalance: { // Aqueous components for (index_t k: m_data->range_aqueous_component()) { const index_t idc = ideq_paq(k); if (idc == no_equation) continue; scalar_t tmp_drdb = 0.0; if (m_data->charge_component(k) != 0.0) tmp_drdb = m_data->charge_component(k); // Secondary species for (index_t j: m_data->range_aqueous()) { if ( not is_aqueous_active(j) or m_data->nu_aqueous(j, k) == 0.0 or m_data->charge_aqueous(j) == 0.0 ) continue; scalar_t tmp_value = m_data->nu_aqueous(j, k)*m_data->charge_aqueous(j); tmp_value *= secondary_molality(j)/component_molality(x, k); tmp_drdb += tmp_value; } jacobian(idp, idc) += component_molality(x,k)*log10*tmp_drdb; } break; } // Fixed activity equation // ======================= case AqueousComponentEquationType::FixedActivity: jacobian(idp, idp) = -1.0/fixed_activity_bc(i); break; // Fixed fugacity equation // ======================= case AqueousComponentEquationType::FixedFugacity: { index_t id_g = m_equations.fixed_activity_species[i]; for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation or m_data->nu_gas(id_g, k) == 0.0) continue; jacobian(idp, ideq_paq(k)) = -m_data->nu_gas(id_g, k)/fixed_fugacity_bc(i); } break; } // end case } // end switch } } void AdimensionalSystem::jacobian_minerals(Vector& x, Matrix& jacobian) { for (index_t m: m_data->range_mineral()) { const index_t idm = ideq_min(m); if (idm == no_equation) continue; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) continue; jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i); } if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m)) jacobian(idm, ideq_electron()) = -m_data->nu_mineral(m, m_data->electron_index()); } } void AdimensionalSystem::jacobian_surface(Vector& x, Matrix& jacobian) { const index_t ids = ideq_surf(); const scalar_t factor = surface_total_concentration(); const scalar_t conc_w = density_water()*saturation_water(x); specmicp_assert(ids != no_equation); scalar_t tmp_s = - free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; tmp_s -= m_data->nb_sorption_sites(s)* m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } jacobian(ids, ids) = conc_w*log10 * tmp_s / factor; // water const index_t idw = ideq_w(); if (idw != no_equation) { const scalar_t rho_w = density_water(); scalar_t tmp_w = - free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; tmp_w -= m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } jacobian(ids, idw) = rho_w * tmp_w / factor; } // component for (index_t k: m_data->range_aqueous_component()) { const index_t idk = ideq_paq(k); if (idk == no_equation) continue; scalar_t tmp_k = - conc_w*diff_surface_weigthed_sum_sorbed(k); jacobian(ids, idk) = tmp_k/factor; } if (ideq_electron() != no_equation) { scalar_t tmp_e = - conc_w*diff_surface_weigthed_sum_sorbed(m_data->electron_index()); jacobian(ids, ideq_electron()) = tmp_e/factor; } // water } void AdimensionalSystem::jacobian_electron(Vector& x, Matrix& jacobian) { const auto ide = ideq_electron(); const auto dofe = m_data->electron_index(); const scalar_t conc_w = density_water()*saturation_water(x); const scalar_t factor = get_options().scaling_electron; // Aqueous components for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp_eip = 0; tmp_eip -= diff_weigthed_sum_aqueous(k, dofe); tmp_eip -= diff_weigthed_sum_sorbed(k, dofe); tmp_eip *= conc_w; tmp_eip -= diff_weigthed_sum_gas(k, dofe); jacobian(ide, ideq_paq(k)) = tmp_eip/factor; } // Minerals for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(ide, ideq_min(m)) = - m_data->nu_mineral(m, dofe)/molar_volume_mineral(m)/factor; } // Water if (ideq_w() != no_equation) { scalar_t tmp_iw = 0; tmp_iw -= weigthed_sum_aqueous(dofe); tmp_iw -= weigthed_sum_sorbed(dofe); tmp_iw *= density_water(); jacobian(ide, ideq_w()) = tmp_iw/factor; } // Surface if (ideq_surf() != no_equation) { scalar_t tmp_s = -conc_w*diff_surface_weigthed_sum_sorbed(dofe); jacobian(ide, ideq_surf()) = tmp_s/factor; } // Electron if (ideq_electron() != no_equation) { scalar_t tmp = 0.0; tmp -= diff_weigthed_sum_aqueous(dofe, dofe); tmp -= diff_weigthed_sum_sorbed(dofe, dofe); tmp*= conc_w; tmp -= diff_weigthed_sum_gas(dofe, dofe); jacobian(ide, ide) = tmp/factor; } } // ========================== // // // // Secondary variables // // // // ========================== // void AdimensionalSystem::set_secondary_variables(const Vector& x) { set_saturation_gas_phase(x); set_pressure_fugacity(x); if (ideq_surf() != no_equation) set_sorbed_concentrations(x); set_secondary_concentration(x); compute_log_gamma(x); } void AdimensionalSystem::set_saturation_gas_phase(const Vector& x) { m_second.saturation_gas = 1 - saturation_water(x) - sum_saturation_minerals(x) - m_inert_volume_fraction; } void AdimensionalSystem::set_pressure_fugacity(const Vector& x) { const auto rt = constants::gas_constant*temperature(); for (index_t k: m_data->range_gas()) { if (not is_active_gas(k)) continue; scalar_t logp = -m_data->logk_gas(k); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_gas(k, i) == 0.0) continue; const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i); logp += m_data->nu_gas(k, i) * log_activity_i; } if (ideq_electron() != no_equation and m_data->is_gas_half_cell_reaction(k)) logp += m_data->nu_gas(k, m_data->electron_index())*log_activity_electron(x); m_second.gas_fugacity(k) = pow10(logp); - const auto pressure = gas_fugacity(k)*gas_total_pressure(); - const auto concentration = saturation_gas_phase()*pressure/rt; + const scalar_t pressure = gas_fugacity(k)*gas_total_pressure(); + const scalar_t concentration = saturation_gas_phase()*pressure/rt; m_second.gas_concentration(k) = concentration; } } void AdimensionalSystem::set_secondary_concentration(const Vector& x) { for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) { m_second.secondary_molalities(j) = 0.0; continue; } scalar_t logconc = -m_data->logk_aqueous(j) - log_gamma_secondary(j); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_aqueous(j, k) == 0) continue; const auto log_activity_k = log_component_molality(x, k) + log_gamma_component(k); logconc += m_data->nu_aqueous(j, k)*log_activity_k; } if (ideq_electron() != no_equation and m_data->is_half_cell_reaction(j)) logconc += m_data->nu_aqueous(j, m_data->electron_index())*log_activity_electron(x); m_second.secondary_molalities(j) = pow10(logconc); } } void AdimensionalSystem::set_sorbed_concentrations(const Vector& x) { for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) { m_second.sorbed_concentrations(s) = 0.0; continue; } scalar_t logconc = -m_data->logk_sorbed(s) + m_data->nb_sorption_sites(s)*(log_free_sorption_site_concentration(x)); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_sorbed(s, k) == 0.0) continue; const auto activity_k = log_component_molality(x, k) + log_gamma_component(k); logconc += m_data->nu_sorbed(s, k)*activity_k; } if (ideq_electron() != no_equation and m_data->is_sorbed_half_cell_reaction(s)) logconc += m_data->nu_sorbed(s, m_data->electron_index())*log_activity_electron(x); m_second.sorbed_concentrations(s) = pow10(logconc); } } void AdimensionalSystem::set_ionic_strength(const Vector& x) { scalar_t ionic = 0; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation or m_data->charge_component(i) == 0) continue; ionic += component_molality(x, i)*std::pow(m_data->charge_component(i),2); } for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j) or m_data->charge_aqueous(j) == 0) continue; ionic += secondary_molality(j)*std::pow(m_data->charge_aqueous(j),2); } ionic_strength() = ionic/2; } void AdimensionalSystem::compute_log_gamma(const Vector& x) { set_ionic_strength(x); const scalar_t sqrti = std::sqrt(ionic_strength()); for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) { log_gamma_component(i) = 0.0; continue; } log_gamma_component(i) = laws::extended_debye_huckel( ionic_strength(), sqrti, m_data->charge_component(i), m_data->a_debye_component(i), m_data->b_debye_component(i) ); } for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) { log_gamma_secondary(j) = 0.0; continue; } log_gamma_secondary(j) = laws::extended_debye_huckel( ionic_strength(), sqrti, m_data->charge_aqueous(j), m_data->a_debye_aqueous(j), m_data->b_debye_aqueous(j) ); } } bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual) { if (not get_options().non_ideality) { set_secondary_variables(x); // we still need to compute secondary species ! // if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) // { // set_secondary_variables(x); // } return true; } not_in_linesearch = true; scalar_t previous_norm = m_second.loggamma.norm(); bool may_have_converged = false; if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) { // Use fixed point iterations for non-ideality for (int i=0; i 1e-6) { WARNING << "Activity coefficient have not converged !" << std::endl << "output can not be trusted\n Difference : " +std::to_string(std::abs(previous_norm - m_second.loggamma.norm())); } } // Set the correct value for the water total saturation if (ideq_w() == no_equation) { xtot(dof_water()) = saturation_water(x); } return AdimensionalSystemSolution(xtot, m_second.secondary_molalities, m_second.loggamma, m_second.ionic_strength, m_second.gas_fugacity, m_second.sorbed_concentrations, m_inert_volume_fraction); } // Water, saturation and density // ============================== scalar_t AdimensionalSystem::density_water() const { return laws::density_water(units::celsius(25.0), length_unit(), mass_unit()); } scalar_t AdimensionalSystem::saturation_water(const Vector& x) const { if (ideq_w() != no_equation) return x(ideq_w()); else return 1.0 - sum_saturation_minerals(x) - m_inert_volume_fraction; } scalar_t AdimensionalSystem::saturation_mineral(const Vector& x, index_t mineral) const { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral()); if (ideq_min(mineral) == no_equation) return 0.0; else return x(ideq_min(mineral)); } scalar_t AdimensionalSystem::sum_saturation_minerals(const Vector& x) const { scalar_t sum_saturations = 0.0; for (index_t mineral: m_data->range_mineral()) { sum_saturations += saturation_mineral(x, mineral); } return sum_saturations; } // Starting guess // ============== void AdimensionalSystem::reasonable_starting_guess(Vector &xtot) { xtot.resize(total_dofs()); xtot(dof_water()) = 0.8; for (index_t i: m_data->range_aqueous_component()) { xtot(dof_component(i)) = -4.0; } if (ideq_surf() != no_equation) xtot(dof_surface()) = std::log10(0.8*surface_total_concentration()); else xtot(dof_surface()) = -HUGE_VAL; if (ideq_electron() != no_equation) xtot(dof_electron()) = -4; else xtot(dof_electron()) = -HUGE_VAL; xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero(); m_second = SecondaryVariables(m_data); } void AdimensionalSystem::reasonable_restarting_guess(Vector& xtot) { std::random_device rd; std::mt19937 gen(rd()); std::uniform_real_distribution<> dis(-2, 2); xtot(dof_water()) = 0.5; for (index_t i: m_data->range_aqueous_component()) { //if (xtot(dof_component(i)) > 0 or xtot(dof_component(i)) < -9) xtot(i) = get_options().restart_concentration + dis(gen); } if (ideq_surf() != no_equation) xtot(dof_surface()) = std::log10(0.8*surface_total_concentration()); else xtot(dof_surface()) = -HUGE_VAL; if (ideq_electron() != no_equation) xtot(dof_electron()) = -4; else xtot(dof_electron()) = -HUGE_VAL; xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero(); m_second = SecondaryVariables(m_data); } } // end namespace specmicp