diff --git a/src/database/data_container.cpp b/src/database/data_container.cpp
index 8bdad54..eb3b667 100644
--- a/src/database/data_container.cpp
+++ b/src/database/data_container.cpp
@@ -1,65 +1,75 @@
/*-------------------------------------------------------
- Module : database
- File : data_container.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#include "data_container.hpp"
namespace specmicp {
namespace database {
//! \brief Return the molar mass (kg/mol) of a mineral
scalar_t DataContainer::molar_mass_mineral(index_t m) {
- return 1e-3*nu_mineral.row(m).dot(molar_mass_basis);
+ return 1e-3*nu_mineral.row(m).dot(_molar_mass_basis);
}
//! \brief Return the molar mass (kg/mol) of a mineral governed by kinetic
scalar_t DataContainer::molar_mass_mineral_kinetic(index_t m) {
- return 1e-3*nu_mineral_kinetic.row(m).dot(molar_mass_basis);
+ return 1e-3*nu_mineral_kinetic.row(m).dot(_molar_mass_basis);
+}
+
+scalar_t DataContainer::molar_volume_mineral(index_t m, units::LengthUnit length_unit)
+{
+ if (_molar_volume_mineral(m) < 0) {
+ throw db_noninitialized_value("molar volume", labels_minerals[m]);}
+ if (length_unit == units::LengthUnit::meter)
+ return 1e-6*_molar_volume_mineral(m);
+ else // cemtimeter
+ return _molar_volume_mineral(m);
}
//! \brief Return the molar volume (m^3/mol) of a mineral
scalar_t DataContainer::molar_volume_mineral(index_t m)
{
if (_molar_volume_mineral(m) < 0) {
throw db_noninitialized_value("molar volume", labels_minerals[m]);}
return 1e-6*_molar_volume_mineral(m);
}
//! \brief Return the molar volume (m^3/mol) of a mineral governed by kinetic
scalar_t DataContainer::molar_volume_mineral_kinetic(index_t m)
{
if (_molar_volume_mineral_kinetic(m) < 0) {
throw db_noninitialized_value("molar volume", labels_minerals_kinetic[m]);}
return 1e-6*_molar_volume_mineral_kinetic(m);
}
} // end namespace database
} // end namespace specmicp
diff --git a/src/database/data_container.hpp b/src/database/data_container.hpp
index 306ca49..302f463 100644
--- a/src/database/data_container.hpp
+++ b/src/database/data_container.hpp
@@ -1,163 +1,170 @@
/*-------------------------------------------------------
- Module : database
- File : data_container.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_DATABASE_DATACONTAINER_HPP
#define SPECMICP_DATABASE_DATACONTAINER_HPP
//! \file data_container.hpp Storage class for thermodynamics database
#include "common_def.hpp"
#include <map>
#include "boost/range/irange.hpp"
+#include "physics/units.hpp"
+
namespace specmicp {
namespace database {
//! \brief Storage class - Contains the database
//!
//! - Should not be accessed directly but through interfaces
//! Correct interfaces are subclasses of DatabaseModule
//! - Should be shared with a smart pointer
struct DataContainer
{
DataContainer():
is_canonical(false)
{}
// Basis
// ------
index_t nb_component; //!< Number of components == size of the basis
vector_labels_t labels_basis; //!< labels of the components
std::map<std::string, index_t> map_labels_basis; //!< map labels <-> id
- data_vector_t molar_mass_basis; //!< molar mass of the basis // (g/mol)
+ data_vector_t _molar_mass_basis; //!< molar mass of the basis // (g/mol)
+ scalar_t molar_mass_basis(index_t component) {return _molar_mass_basis(component);}
scalar_t molar_mass_basis_si(index_t i) {return 1e-3*molar_mass_basis(i);}
Eigen::Matrix<scalar_t, Eigen::Dynamic, 3> param_aq; //!< Aqueous parameters (both for basis and secondary species)
//! Return the charge of a component
scalar_t& charge_component(index_t i) {return param_aq(i, 0);}
//! Return the 'a_i' Debye-Huckel parameter for a component
scalar_t& a_debye_component(index_t i) {return param_aq(i, 1);}
//! Return the 'b_i' extended Debye-Huckel parameter for a component
scalar_t& b_debye_component(index_t i) {return param_aq(i, 2);}
// Secondary aqueous species
// -------------------------
index_t nb_aqueous; //!< Number of aqueous species (not taking into acount the basis)
vector_labels_t labels_aqueous; //!< labels of the aqueous species
reaction_mat_t nu_aqueous; //!< Stoechiometric coefficient for aqueous species
logK_vector_t logk_aqueous; //!< LogK for aqueous species
//! Return the charge of a secondary aqueous species
scalar_t& charge_aqueous(index_t j) {return param_aq(j+nb_component, 0);}
//! Return the 'a_i' Debye-Huckel parameter for a secondary aqueous species
scalar_t& a_debye_aqueous(index_t j) {return param_aq(j+nb_component, 1);}
//! Return the 'b_i' extended Debye-Huckel parameter for a secondary aqueous species
scalar_t& b_debye_aqueous(index_t j) {return param_aq(j+nb_component, 2);}
// Minerals
// ---------
index_t nb_mineral; //!< Number of minerals (used in computation)
vector_labels_t labels_minerals; //!< labels of the minerals
reaction_mat_t nu_mineral; //!< Stoichiometric coefficient for minerals
logK_vector_t logk_mineral; //!< LogK for minerals
data_vector_t _molar_volume_mineral; //!< molar volume of mineral (cm3/mol)
list_stability_t _stability_mineral; //!< Stability of a mineral
MineralStabilityClass stability_mineral(index_t m) const {return _stability_mineral[m];}
MineralStabilityClass& stability_mineral(index_t m) {return _stability_mineral[m];}
index_t nb_mineral_kinetic; //!< Number of minerals (governed by kinetics)
vector_labels_t labels_minerals_kinetic; //!< labels of the minerals (governed by kinetics)
reaction_mat_t nu_mineral_kinetic; //!< Stoichiometric coefficient for minerals (governed by kinetics)
logK_vector_t logk_mineral_kinetic; //!< LogK for minerals (governed by kinetics)
data_vector_t _molar_volume_mineral_kinetic; //!< molar volume of mineral (cm3/mol)
// Gas
// ----
index_t nb_gas; //!< Number of gas
vector_labels_t labels_gas; //!< labels of the gas phase
reaction_mat_t nu_gas; //!< Stoichiometric coefficient for the gas
logK_vector_t logk_gas; //!< log(K) gas
// Ranges
// ------
range_t range_component() { return boost::irange((index_t) 0, nb_component);}
range_t range_aqueous_component() { return boost::irange((index_t) 1, nb_component);}
range_t range_aqueous() { return boost::irange((index_t) 0, nb_aqueous);}
range_t range_mineral() { return boost::irange((index_t) 0, nb_mineral);}
range_t range_mineral_kinetic() { return boost::irange((index_t) 0, nb_mineral_kinetic);}
range_t range_gas() { return boost::irange((index_t) 0, nb_gas);}
// Status
// ------
bool is_canonical; //!< true if the database is in canonical form
// Computed variables
// ------------------
// The following functions are provided for convenience
// They provide variables computed from the main data
//! \brief Return the molar mass (kg/mol) of a mineral
scalar_t molar_mass_mineral(index_t m);
//! \brief Return the molar mass (kg/mol) of a mineral governed by kinetic
scalar_t molar_mass_mineral_kinetic(index_t m);
//! \brief Return the molar volume (m^3/mol) of a mineral
scalar_t molar_volume_mineral(index_t m);
//! \brief Return the molar volume (m^3/mol) of a mineral governed by kinetic
scalar_t molar_volume_mineral_kinetic(index_t m);
+ //! Return the molar volume of a mineral in mol/(length_unit)^3
+ scalar_t molar_volume_mineral(index_t m, units::LengthUnit length_unit);
+ //! Return the molar mass of 'component' in mass_unit/mol
+ scalar_t molar_mass_basis(index_t component, units::MassUnit mass_unit);
};
} // end namespace database
} // end namespace specmicp
#include <memory>
namespace specmicp {
namespace database {
//! \brief type of the database
using RawDatabasePtr = std::shared_ptr<DataContainer>;
} // end namespace database
} // end namespace specmicp
#endif // SPECMICP_DATABASE_DATACONTAINER_HPP
diff --git a/src/database/reader.cpp b/src/database/reader.cpp
index 748057c..c71b2c2 100644
--- a/src/database/reader.cpp
+++ b/src/database/reader.cpp
@@ -1,497 +1,497 @@
/*-------------------------------------------------------
- Module : database
- File : reader.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#include "reader.hpp"
#include "json/json.h"
#include <fstream>
#include <stdexcept>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/trim_all.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <locale>
#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_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_FLAG_KINETIC "flag_kinetic"
#define INDB_OPEN_DELIMITER_CHARGE '['
#define INDB_CLOSE_DELIMITER_CHARGE ']'
namespace specmicp {
namespace database {
void DataReader::parse()
{
std::ifstream datafile(m_filepath, std::ifstream::binary);
Json::Reader reader;
if (not datafile.is_open())
{
ERROR << "Database not found : " << m_filepath;
std::string message("Database not found : " + m_filepath);
throw std::invalid_argument(message);
}
DEBUG << "Parsing database";
bool parsingSuccessful = reader.parse( datafile, m_root);
if ( !parsingSuccessful )
{
// report to the user the failure and their locations in the document.
std::string message("Failed to parse configuration\n" + reader.getFormatedErrorMessages());
throw std::runtime_error(message);
}
parse_size_db();
parse_basis();
parse_aqueous();
parse_minerals();
parse_gas();
datafile.close();
}
void DataReader::check_database()
{
// TODO
}
void DataReader::parse_size_db()
{
data->nb_component = m_root[INDB_SECTION_BASIS].size();
data->nb_aqueous = m_root[INDB_SECTION_AQUEOUS].size();
data->nb_mineral = m_root[INDB_SECTION_MINERALS].size();
data->nb_gas = m_root[INDB_SECTION_GAS].size();
data->param_aq = Eigen::MatrixXd::Zero(data->nb_component+data->nb_aqueous, 3);
}
void DataReader::parse_basis()
{
Json::Value& basis = get_mandatory_section(m_root, INDB_SECTION_BASIS);
DEBUG << "Size basis : " << basis.size();
int size_basis = basis.size();
data->labels_basis.reserve(size_basis);
- data->molar_mass_basis.resize(size_basis);
+ data->_molar_mass_basis.resize(size_basis);
for (int id=0; id<size_basis; ++id)
{
Json::Value& species = basis[id];
check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL);
std::string label = species[INDB_ATTRIBUTE_LABEL].asString();
SPAM << "Parsing component : " << label;
data->map_labels_basis[label] = id;
data->labels_basis.push_back(label);
data->param_aq(id, 0) = charge_from_label(label);
if (species.isMember(INDB_ATTRIBUTE_ACTIVITY))
{
data->param_aq(id, 1) = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_A].asDouble();
data->param_aq(id, 2) = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_B].asDouble();
}
check_for_mandatory_value(species, INDB_ATTRIBUTE_MOLARMASS);
- data->molar_mass_basis(id) = species[INDB_ATTRIBUTE_MOLARMASS].asDouble();
+ data->_molar_mass_basis(id) = species[INDB_ATTRIBUTE_MOLARMASS].asDouble();
}
}
void DataReader::parse_aqueous()
{
DEBUG << "Parse aqueous species";
Json::Value& aqueous = get_mandatory_section(m_root, INDB_SECTION_AQUEOUS);
data->labels_aqueous.reserve(data->nb_aqueous);
data->nu_aqueous = reaction_mat_t::Zero(data->nb_aqueous, data->nb_component+data->nb_aqueous);
data->logk_aqueous = logK_vector_t(data->nb_aqueous);
for (int id=0; id<data->nb_aqueous; ++id)
{
Json::Value& species = aqueous[id];
check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL);
std::string label = species[INDB_ATTRIBUTE_LABEL].asString();
data->labels_aqueous.push_back(label);
const int idparam = id + data->nb_component;
data->param_aq(idparam, 0) = charge_from_label(label);
if (species.isMember(INDB_ATTRIBUTE_ACTIVITY))
{
data->param_aq(idparam, 1) = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_A].asDouble();
data->param_aq(idparam, 2) = species[INDB_ATTRIBUTE_ACTIVITY][INDB_ATTRIBUTE_ACTIVITY_B].asDouble();
}
// equation
check_for_mandatory_value(species, INDB_ATTRIBUTE_COMPOSITION);
std::map<std::string, scalar_t> compo;
parse_equation(species[INDB_ATTRIBUTE_COMPOSITION].asString(), compo);
double test_charge=0;
for (auto it=compo.begin(); it != compo.end(); ++it)
{
std::map<std::string, index_t>::const_iterator search = data->map_labels_basis.find(it->first);
if(search != data->map_labels_basis.end()) {
data->nu_aqueous(id, search->second) = it->second;
test_charge += it->second*data->param_aq(search->second, 0);
}
else {
vector_labels_t::const_iterator searchaq = std::find(data->labels_aqueous.begin(),
data->labels_aqueous.end(),
it->first);
if (searchaq != data->labels_aqueous.end())
{
const int id_aq = (searchaq - data->labels_aqueous.begin()) + data->nb_component;
data->nu_aqueous(id, id_aq) = it->second;
test_charge += it->second*data->param_aq(id_aq, 0);
}
else
{
throw db_invalid_syntax("Unknown species : '" + it->first + "' while parsing equations for " + label +".");
}
}
}
if (test_charge != data->param_aq(idparam, 0))
{
throw db_invalid_data("Total charge is not zero for aqueous species : '"+label+
"' - Charge-decomposition - charge species : "+
std::to_string(test_charge - data->param_aq(id, 0) )+".");
}
check_for_mandatory_value(species, INDB_ATTRIBUTE_LOGK);
data->logk_aqueous(id) = species[INDB_ATTRIBUTE_LOGK].asDouble();
}
}
void DataReader::parse_minerals()
{
DEBUG << "Parse minerals";
Json::Value& minerals = get_mandatory_section(m_root, INDB_SECTION_MINERALS);
// find kinetic minerals
int nb_kin = 0;
for (int id=0; id<data->nb_mineral; ++id)
{
bool is_kin = false;
if (minerals[id].isMember(INDB_ATTRIBUTE_FLAG_KINETIC))
{
is_kin = minerals[id][INDB_ATTRIBUTE_FLAG_KINETIC].asBool();
}
if (is_kin) ++nb_kin;
}
data->nb_mineral = data->nb_mineral - nb_kin;
data->labels_minerals.reserve(data->nb_mineral);
data->nu_mineral = reaction_mat_t::Zero(data->nb_mineral, data->nb_component+data->nb_aqueous);
data->logk_mineral = logK_vector_t(data->nb_mineral);
data->_molar_volume_mineral.resize(data->nb_mineral);
data->_stability_mineral = list_stability_t(data->nb_mineral, MineralStabilityClass::equilibrium);
data->nb_mineral_kinetic = nb_kin;
data->labels_minerals_kinetic.reserve(data->nb_mineral_kinetic);
data->nu_mineral_kinetic = reaction_mat_t::Zero(data->nb_mineral_kinetic, data->nb_component+data->nb_aqueous);
data->logk_mineral_kinetic = logK_vector_t(data->nb_mineral_kinetic);
data->_molar_volume_mineral_kinetic.resize(data->nb_mineral);
// id for each class of minerals
int id_in_eq=0;
int id_in_kin=0;
for (int id=0; id<data->nb_mineral+data->nb_mineral_kinetic; ++id)
{
Json::Value& species = minerals[id];
//check if material is at equilibrium or governed by kinetics
bool is_kin = false;
if (species.isMember(INDB_ATTRIBUTE_FLAG_KINETIC))
{
is_kin = species[INDB_ATTRIBUTE_FLAG_KINETIC].asBool();
}
check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL);
const std::string label = species[INDB_ATTRIBUTE_LABEL].asString();
if (is_kin) data->labels_minerals_kinetic.push_back(label);
else data->labels_minerals.push_back(label);
// equation
// --------
check_for_mandatory_value(species, INDB_ATTRIBUTE_COMPOSITION);
std::map<std::string, scalar_t> compo;
parse_equation(species[INDB_ATTRIBUTE_COMPOSITION].asString(), compo);
double test_charge = 0;
for (auto it=compo.begin(); it != compo.end(); ++it)
{
std::map<std::string, index_t>::const_iterator search = data->map_labels_basis.find(it->first);
if(search != data->map_labels_basis.end()) { // this is a component
if (is_kin) data->nu_mineral_kinetic(id_in_kin, search->second) = it->second;
else data->nu_mineral(id_in_eq, search->second) = it->second;
test_charge += it->second * data->param_aq(search->second, 0);
}
else { // this is an aqueous species
vector_labels_t::const_iterator searchaq = std::find(data->labels_aqueous.begin(),
data->labels_aqueous.end(),
it->first);
if (searchaq != data->labels_aqueous.end())
{
int id_aq = searchaq-data->labels_aqueous.begin()+data->nb_component;
if (is_kin) data->nu_mineral_kinetic(id_in_kin,id_aq) = it->second;
else data->nu_mineral(id_in_eq, id_aq) = it->second;
test_charge += it->second * data->param_aq(id_aq, 0);
}
else // or something we don't know
{
throw db_invalid_syntax("Unknown species : '" + it->first + "' while parsing equations for " + label +".");
}
}
}
if (test_charge != 0)
{
throw db_invalid_data("Total charge is not zero for mineral : '"+label+
"' - total charge : "+std::to_string(test_charge)+".");
}
// log(K)
// ------
check_for_mandatory_value(species, INDB_ATTRIBUTE_LOGK);
if (is_kin) data->logk_mineral_kinetic(id_in_kin) = species[INDB_ATTRIBUTE_LOGK].asDouble();
else data->logk_mineral(id_in_eq) = species[INDB_ATTRIBUTE_LOGK].asDouble();
// Molar volume
// -------------
if (species.isMember(INDB_ATTRIBUTE_MOLARVOLUME))
{
if (is_kin) data->_molar_volume_mineral_kinetic(id_in_kin) = species[INDB_ATTRIBUTE_MOLARVOLUME].asDouble();
else data->_molar_volume_mineral(id_in_eq) = species[INDB_ATTRIBUTE_MOLARVOLUME].asDouble();
}
else
{
if (is_kin) data->_molar_volume_mineral_kinetic(id_in_kin) = -1;
else data->_molar_volume_mineral(id_in_eq) = -1;
}
// update indices
if (is_kin) {++id_in_kin;}
else {++id_in_eq;}
}
assert(id_in_kin == data->nb_mineral_kinetic);
assert(id_in_eq == data->nb_mineral);
}
void DataReader::parse_gas()
{
DEBUG << "Parse gas";
Json::Value& gas = get_mandatory_section(m_root, INDB_SECTION_GAS);
data->labels_gas.reserve(data->nb_gas);
data->nu_gas = reaction_mat_t::Zero(data->nb_gas, data->nb_component+data->nb_aqueous);
data->logk_gas = logK_vector_t(data->nb_gas);
for (int id=0; id<data->nb_gas; ++id)
{
Json::Value& species = gas[id];
check_for_mandatory_value(species, INDB_ATTRIBUTE_LABEL);
std::string label = species[INDB_ATTRIBUTE_LABEL].asString();
data->labels_gas.push_back(label);
// equation
check_for_mandatory_value(species, INDB_ATTRIBUTE_COMPOSITION);
std::map<std::string, scalar_t> compo;
parse_equation(species[INDB_ATTRIBUTE_COMPOSITION].asString(), compo);
double test_charge=0;
for (auto it=compo.begin(); it != compo.end(); ++it)
{
std::map<std::string, index_t>::const_iterator search = data->map_labels_basis.find(it->first);
if(search != data->map_labels_basis.end()) {
data->nu_gas(id, search->second) = it->second;
test_charge += it->second*data->param_aq(search->second, 0);
}
else {
vector_labels_t::const_iterator searchaq = std::find(data->labels_aqueous.begin(),
data->labels_aqueous.end(),
it->first);
if (searchaq != data->labels_aqueous.end())
{
const int id_aq = (searchaq - data->labels_aqueous.begin()) + data->nb_component;
data->nu_gas(id, id_aq) = it->second;
test_charge += it->second*data->param_aq(id_aq, 0);
}
else
{
throw db_invalid_syntax("Unknown species : '" + it->first + "' while parsing equations for " + label +".");
}
}
}
if (test_charge != 0)
{
throw db_invalid_data("Total charge is not zero for gas '"+label+
"' : " + std::to_string(test_charge)+".");
}
check_for_mandatory_value(species, INDB_ATTRIBUTE_LOGK);
data->logk_gas(id) = species[INDB_ATTRIBUTE_LOGK].asDouble();
}
}
void parse_equation(const std::string& equation, std::map<std::string, scalar_t> &compo)
{
std::vector<std::string> 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;
}
} // end namespace specmicp
} // end namespace chemy
#undef INDB_SECTION_BASIS
#undef INDB_SECTION_AQUEOUS
#undef INDB_SECTION_MINERALS
#undef INDB_SECTION_GAS
#undef INDB_ATTRIBUTE_NAME
#undef INDB_ATTRIBUTE_LABEL
#undef INDB_ATTRIBUTE_ACTIVITY
#undef INDB_ATTRIBUTE_ACTIVITY_A
#undef INDB_ATTRIBUTE_ACTIVITY_B
#undef INDB_ATTRIBUTE_MOLARMASS
#undef INDB_ATTRIBUTE_COMPOSITION
#undef INDB_ATTRIBUTE_LOGK
#undef INDB_ATTRIBUTE_FLAG_KINETIC
#undef INDB_OPEN_DELIMITER_CHARGE
#undef INDB_CLOSE_DELIMITER_CHARGE
diff --git a/src/database/selector.cpp b/src/database/selector.cpp
index 29b74f3..5c1d73a 100644
--- a/src/database/selector.cpp
+++ b/src/database/selector.cpp
@@ -1,280 +1,280 @@
/*-------------------------------------------------------
- Module : database
- File : selector.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#include "selector.hpp"
namespace specmicp {
namespace database {
void DatabaseSelector::remove_component(const std::vector<index_t>& id_component_to_remove)
{
analyse_component(id_component_to_remove);
select_aqueous(id_component_to_remove);
select_minerals(id_component_to_remove);
select_minerals_kinetic(id_component_to_remove);
if (data->nb_gas != 0) select_gas(id_component_to_remove);
select_components();
}
void DatabaseSelector::analyse_component(const std::vector<index_t>& id_component_to_remove)
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
m_is_component_to_remove[*it] = true;
}
m_nb_component_to_keep = data->nb_component - id_component_to_remove.size();
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->param_aq.row(new_i) = data->param_aq.row(i);
- data->molar_mass_basis(new_i) = data->molar_mass_basis(i);
+ data->_molar_mass_basis(new_i) = data->molar_mass_basis(i);
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
- data->molar_mass_basis.conservativeResize(nb_component_to_keep());
+ data->_molar_mass_basis.conservativeResize(nb_component_to_keep());
}
void DatabaseSelector::select_aqueous(const std::vector<index_t>& id_component_to_remove)
{
// first we select which aqueous species we should remove
std::vector<bool> aqueous_to_remove(data->nb_aqueous, false);
for (index_t j: data->range_aqueous())
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_aqueous(j, *it) != 0)
{
aqueous_to_remove[j] = true;
break;
}
}
}
// then we remove them, in two steps
const uindex_t nb_aq_to_keep = std::count(aqueous_to_remove.begin(), aqueous_to_remove.end(), false);
uindex_t new_j = 0;
// 1) first we copy data in the beginning of the arrays
for (index_t j: data->range_aqueous())
{
if (aqueous_to_remove[j]) continue;
data->logk_aqueous(new_j) = data->logk_aqueous(j);
data->param_aq.row(nb_component_to_keep()+new_j) = data->param_aq.row(data->nb_component+j);
data->labels_aqueous[new_j] = data->labels_aqueous[j];
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->nu_aqueous(new_j, new_i) = data->nu_aqueous(j, i);
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
++new_j;
}
specmicp_assert(new_j == nb_aq_to_keep);
// 2) then we resize the arrays
data->nu_aqueous.conservativeResize(nb_aq_to_keep, nb_component_to_keep());
data->logk_aqueous.conservativeResize(nb_aq_to_keep);
data->labels_aqueous.resize(nb_aq_to_keep);
data->param_aq.conservativeResize(nb_component_to_keep()+nb_aq_to_keep, Eigen::NoChange);
data->nb_aqueous = nb_aq_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_minerals(const std::vector<index_t>& id_component_to_remove)
{
// first we select which minerals we should remove
std::vector<bool> minerals_to_remove(data->nb_mineral, false);
for (index_t m: data->range_mineral())
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_mineral(m, *it) != 0)
{
minerals_to_remove[m] = true;
break;
}
}
}
// then we remove them, in two steps
const uindex_t nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
uindex_t new_m = 0;
// 1) first we copy data in the beginning of the arrays
for (index_t m: data->range_mineral())
{
if (minerals_to_remove[m]) continue;
data->logk_mineral(new_m) = data->logk_mineral(m);
data->labels_minerals[new_m] = data->labels_minerals[m];
data->_molar_volume_mineral(new_m) = data->_molar_volume_mineral(m);
data->stability_mineral(new_m) = data->stability_mineral(m);
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->nu_mineral(new_m, new_i) = data->nu_mineral(m, i);
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
++new_m;
}
specmicp_assert(new_m == nb_min_to_keep);
// 2) then we resize the arrays
data->nu_mineral.conservativeResize(nb_min_to_keep, nb_component_to_keep());
data->logk_mineral.conservativeResize(nb_min_to_keep);
data->labels_minerals.resize(nb_min_to_keep);
data->_molar_volume_mineral.conservativeResize(nb_min_to_keep);
data->_stability_mineral.resize(new_m);
data->nb_mineral = nb_min_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_minerals_kinetic(const std::vector<index_t>& id_component_to_remove)
{
// first we select which minerals we should remove
std::vector<bool> minerals_to_remove(data->nb_mineral_kinetic, false);
for (index_t m: data->range_mineral_kinetic())
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_mineral_kinetic(m, *it) != 0)
{
minerals_to_remove[m] = true;
break;
}
}
}
// then we remove them, in two steps
const uindex_t nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
uindex_t new_m = 0;
// 1) first we copy data in the beginning of the arrays
for (index_t m: data->range_mineral_kinetic())
{
if (minerals_to_remove[m]) continue;
data->logk_mineral_kinetic(new_m) = data->logk_mineral_kinetic(m);
data->labels_minerals_kinetic[new_m] = data->labels_minerals_kinetic[m];
data->_molar_volume_mineral_kinetic(new_m) = data->_molar_volume_mineral_kinetic(m);
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->nu_mineral_kinetic(new_m, new_i) = data->nu_mineral_kinetic(m, i);
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
++new_m;
}
specmicp_assert(new_m == nb_min_to_keep);
// 2) then we resize the arrays
data->nu_mineral_kinetic.conservativeResize(nb_min_to_keep, nb_component_to_keep());
data->logk_mineral_kinetic.conservativeResize(nb_min_to_keep);
data->labels_minerals_kinetic.resize(nb_min_to_keep);
data->_molar_volume_mineral_kinetic.conservativeResize(nb_min_to_keep);
data->nb_mineral_kinetic = nb_min_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_gas(const std::vector<index_t>& id_component_to_remove)
{
// first we select which aqueous species we should remove
std::vector<bool> gas_to_remove(data->nb_gas, false);
for (index_t j: data->range_gas())
{
for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
{
if (data->nu_gas(j, *it) != 0)
{
gas_to_remove[j] = true;
break;
}
}
}
// then we remove them, in two steps
const uindex_t nb_gas_to_keep = std::count(gas_to_remove.begin(), gas_to_remove.end(), false);
uindex_t new_j = 0;
// 1) first we copy data in the beginning of the arrays
for (index_t j: data->range_gas())
{
if (gas_to_remove[j]) continue;
data->logk_gas(new_j) = data->logk_gas(j);
data->labels_gas[new_j] = data->labels_gas[j];
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->nu_gas(new_j, new_i) = data->nu_gas(j, i);
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
++new_j;
}
specmicp_assert(new_j == nb_gas_to_keep);
// 2) then we resize the arrays
data->nu_gas.conservativeResize(nb_gas_to_keep, nb_component_to_keep());
data->logk_gas.conservativeResize(nb_gas_to_keep);
data->labels_gas.resize(nb_gas_to_keep);
data->nb_gas = nb_gas_to_keep;
// By doing in this order, we avoid bulk copy of arrays
}
void DatabaseSelector::select_components()
{
uindex_t new_i = 0;
for (index_t i: data->range_component())
{
if (is_component_to_remove(i)) continue;
data->labels_basis[new_i] = data->labels_basis[i];
++new_i;
}
specmicp_assert(new_i == nb_component_to_keep());
data->labels_basis.resize(nb_component_to_keep());
data->nb_component = nb_component_to_keep();
}
void DatabaseSelector::keep_only_component(const std::vector<index_t>& id_to_keep)
{
// First build the list of components to remove
std::vector<index_t> id_to_remove;
id_to_remove.reserve(data->nb_component - id_to_keep.size());
for (index_t id: data->range_component())
{
auto search = std::find(id_to_keep.begin(), id_to_keep.end(), id);
if (search == id_to_keep.end())
id_to_remove.push_back(id);
}
// Then remove them
remove_component(id_to_remove);
}
} // end namespace database
} // end namespace specmicp
diff --git a/src/database/switch_basis.cpp b/src/database/switch_basis.cpp
index f712524..b040a9d 100644
--- a/src/database/switch_basis.cpp
+++ b/src/database/switch_basis.cpp
@@ -1,112 +1,112 @@
/*-------------------------------------------------------
- Module : database
- File : switch_basis.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#include "switch_basis.hpp"
#include <Eigen/Eigen>
#include "canonicalizer.hpp"
namespace specmicp {
namespace database {
void BasisSwitcher::change_basis(std::vector<index_t>& new_basis)
{
if (not is_database_canonical())
{
canonicalize_database(data);
}
Eigen::MatrixXd beta = Eigen::MatrixXd::Zero(data->nb_component, data->nb_component);
Eigen::VectorXd kswap(data->nb_component);
// swap coefficients, build beta matrix and kswap
for (index_t i: data->range_component())
{
if (new_basis[i] >= data->nb_component)
{
index_t j = new_basis[i] - data->nb_component;
beta.row(i) = data->nu_aqueous.row(j);
data->nu_aqueous.row(j) = Vector::Zero(data->nb_component);
data->nu_aqueous(j, i) = 1.0;
kswap(i) = data->logk_aqueous(j);
data->logk_aqueous(j) = 0.0;
}
else
{
beta(i, i) = 1.0;
kswap(i) = 0.0;
}
}
Matrix betainv = beta.inverse();
- data->molar_mass_basis = beta*data->molar_mass_basis;
+ data->_molar_mass_basis = beta*data->_molar_mass_basis;
data->nu_aqueous = data->nu_aqueous * betainv;
data->logk_aqueous -= data->nu_aqueous * kswap;
swap_aq_param(new_basis);
swap_labels(new_basis);
// now do the minerals
data->nu_mineral = data->nu_mineral*betainv;
data->logk_mineral -= data->nu_mineral * kswap;
data->nu_mineral_kinetic = data->nu_mineral_kinetic*betainv;
data->logk_mineral_kinetic -= data->nu_mineral_kinetic * kswap;
// and the gases
if (data->nb_gas != 0)
{
data->nu_gas = data->nu_gas*betainv;
data->logk_gas -= data->nu_gas * kswap;
}
}
void BasisSwitcher::swap_aq_param(std::vector<index_t> new_basis)
{
for (index_t i: data->range_component())
{
if(new_basis[i] >= data->nb_component)
{
const index_t j = new_basis[i];
data->param_aq.row(i).swap(data->param_aq.row(j));
}
}
}
void BasisSwitcher::swap_labels(std::vector<index_t> new_basis)
{
for (index_t i: data->range_component())
{
if (new_basis[i] >= data->nb_component)
{
std::swap(data->labels_basis[i], data->labels_aqueous[new_basis[i]-data->nb_component]);
}
}
}
} // end namespace database
} // end namespace specmicp