Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/equilibrium_data.cpp
b53f240bd6a4reactivetransport

Last Change
Raw File

equilibrium_data.cpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : equilibrium_data.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "equilibrium_data.hpp"
#include "database/module.hpp"
namespace specmicp {
//! return log(IAP) for mineral m
double EquilibriumState::logIAP(int m)
{
double logiap = 0.0;
for (int i=1; i<m_data->nb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral(m, i);
return logiap;
}
double EquilibriumState::logIAP(const std::string& label_mineral)
{
database::DatabaseModule namer(m_data);
int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals);
return logIAP(idm);
}
double EquilibriumState::saturation_index(const std::string& label_mineral)
{
database::DatabaseModule namer(m_data);
int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals);
return logIAP(idm) - m_data->logk_mineral(idm);
}
//! return log(IAP) for mineral m
double EquilibriumState::logIAP_kinetic(int m)
{
double logiap = 0.0;
for (int i=1; i<m_data->nb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral_kinetic(m, i);
return logiap;
}
double EquilibriumState::logIAP_kinetic(const std::string& label_mineral)
{
database::DatabaseModule namer(m_data);
int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals_kinetic);
return logIAP_kinetic(idm);
}
double EquilibriumState::saturation_index_kinetic(const std::string& label_mineral)
{
database::DatabaseModule namer(m_data);
int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals);
return logIAP_kinetic(idm) - m_data->logk_mineral_kinetic(idm);
}
//! Return total concentration of aqueous species for component i
double EquilibriumState::total_aqueous_concentration_component(int i)
{
double totconc = POW10(m_main_variables(i));
for (int j=0; j<m_data->nb_aqueous; ++j)
{
if (m_data->nu_aqueous(j,i) != 0)
totconc += m_data->nu_aqueous(j,i)*m_concentration_secondary(j);
}
return totconc;
}
void EquilibriumState::total_aqueous_concentrations(Eigen::VectorXd& totconc)
{
assert(totconc.rows() == m_data->nb_component);
totconc(0) = m_main_variables(0);
for (int i=1; i<m_data->nb_component; ++i) totconc(i) = total_aqueous_concentration_component(i);
}
double EquilibriumState::pH()
{
// search the information about "H[+]"
database::DatabaseModule namer(m_data);
double id = namer.component_label_to_id("H[+]");
if (id != database::DatabaseModule::no_species) return -m_main_variables(id);
else
{
id = namer.aqueous_label_to_id("H[+]");
assert(id != database::DatabaseModule::no_species);
return -std::log10(m_concentration_secondary(id));
}
}
double EquilibriumState::mass_mineral(int m)
{
return m_data->molar_mass_mineral(m)*m_main_variables(m_data->nb_component+m);
}
} // end namespace specmicp
c4science · Help