Homec4science
Diffusion SpecMiCP / ReactMiCP (specmicp_bin) Browse

specmicp/src/specmicp/equilibrium_data.hpp
f0bae6f199e0specmicp_bin

Last Change
Raw File

equilibrium_data.hpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : equilibrium_data.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_SPECMICP_EQUILIBRIUMDATA_HPP
#define SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP
//! \file equilibrium_data.hpp store results from equilibrium computation
#include "database/database.hpp"
#include "simulation_box.hpp"
#define POW10(m) std::pow(10.0, m)
namespace specmicp {
//! \brief Store the equilibrium state of a simulation
//!
//! This class is used to store the output of an equilibrium solver
class EquilibriumState: public SimulationBox
{
public:
//! \brief Blank state
//!
//! Must be initialized after, should not be used except by internals
EquilibriumState():
m_has_been_properly_initialised(false)
{}
//! \brief Copy constructor
//!
EquilibriumState(const EquilibriumState& other):
m_has_been_properly_initialised(other.m_has_been_properly_initialised),
m_main_variables(other.m_main_variables),
m_concentration_secondary(other.m_concentration_secondary),
m_ionic_strength(other.m_ionic_strength),
m_loggamma(other.m_loggamma),
m_data(other.m_data)
{}
//! \brief Create an equilibrium state
//!
//! This constructor is to be used for an aqueous computation
template <class Derived1, class Derived2, class Derived3>
EquilibriumState(const Eigen::MatrixBase<Derived1>& main_vars,
const Eigen::MatrixBase<Derived2>& secondary,
const Eigen::MatrixBase<Derived3>& loggamma,
double ionic_strength,
std::shared_ptr<database::DataContainer> data):
m_has_been_properly_initialised(true),
m_main_variables(main_vars),
m_concentration_secondary(secondary),
m_ionic_strength(ionic_strength),
m_loggamma(loggamma),
m_data(data)
{}
//! \brief Create an equilibrium state
//!
//! This constructor is to be used for an aqueous computation
EquilibriumState(const Eigen::VectorXd& main_vars,
const Eigen::VectorXd& secondary,
const Eigen::VectorXd& loggamma,
double ionic_strength,
std::shared_ptr<database::DataContainer> data):
m_has_been_properly_initialised(true),
m_main_variables(main_vars),
m_concentration_secondary(secondary),
m_ionic_strength(ionic_strength),
m_loggamma(loggamma),
m_data(data)
{}
//! \brief Create an equilibrium state
//!
//! This constructor is to be used for a computation with gas phases
EquilibriumState(const Eigen::VectorXd& main_vars,
const Eigen::VectorXd& secondary,
const Eigen::VectorXd& loggamma,
double ionic_strength,
const Eigen::VectorXd& pressures,
double volume,
std::shared_ptr<database::DataContainer> data):
SimulationBox(volume),
m_has_been_properly_initialised(true),
m_main_variables(main_vars),
m_concentration_secondary(secondary),
m_pressure_gas(pressures),
m_ionic_strength(ionic_strength),
m_loggamma(loggamma),
m_data(data)
{}
//! \brief Return true if it contains actual data
bool is_valid() {return m_has_been_properly_initialised;}
// Main variables
// ==============
//! \brief Return the mass of free water in the system
double mass_water() const {return m_main_variables(0);}
//! \brief Return the molality of component i
double molality_component(int i) const {return POW10(m_main_variables(i));}
//! \brief Return the number of moles of the mineral at equilibrium 'm'
double moles_mineral(int m) const {
assert(m >= 0 and m <m_data->nb_mineral);
return m_main_variables(m_data->nb_component+m);}
//! \brief Return the mass of a mineral (kg)
double mass_mineral(int m) const;
//! \brief Return the volume of a mineral (V)
double volume_mineral(int m) const {
return m_data->molar_volume_mineral(m)*moles_mineral(m);
}
//! Return the volume of all minerals
double volume_minerals() const;
const Eigen::VectorXd& main_variables() const {return m_main_variables;}
// molality
// ========
//! \brief return the pH
double pH() const;
//! \brief Return the molality of secondary species j
double molality_secondary(int j) const {return m_concentration_secondary(j);}
//! \brief Return a const reference to the vector of molality of the secondary array
const Eigen::VectorXd& molality_secondary() const {return m_concentration_secondary;}
//! \brief Return total concentration of aqueous species for component i
double total_aqueous_concentration_component(int i) const;
//! \brief Compute total aqeuous concentration for each component
void total_aqueous_concentrations(Eigen::VectorXd& totconc) const;
// Activity
// ========
double ionic_strength() const {return m_ionic_strength;}
//! \brief Return a const reference to the activity coefficient vector
const Eigen::VectorXd& activity_coefficient() const {return m_loggamma;}
//! \brief Return activity coefficient of component i
double activity_coefficient_component(int i) const {
assert(i >= 1);
return POW10(m_loggamma(i));
}
//! \brief Return logarithm of activity coefficient of component i
double loggamma_component(int i) const {
assert(i >= 1 and i < m_data->nb_component);
return m_loggamma(i);
}
//! \brief Return activity coefficient of secondary species j
double activity_coefficient_secondary(int j) const {
assert(j >= 0 and j < m_data->nb_aqueous);
return POW10(m_loggamma(j+m_data->nb_component));
}
//! \brief Return the log of the activity coefficient of secondary species j
double loggamma_secondary(int j) const {
assert(j >= 0 and j < m_data->nb_aqueous);
return m_loggamma(j+m_data->nb_component);
}
//! \brief Return activity of component i (i > 1)
double activity_component(int i) const {
assert(i >= 1);
return POW10(m_loggamma(i) + m_main_variables(i));
}
//! \brief Return activity of secondary species j
double activity_secondary(int j) const {
assert(j >= 0 and j < m_data->nb_aqueous);
return POW10(m_loggamma(j+m_data->nb_component))*m_concentration_secondary(j);
}
// Minerals
// ========
//! \brief return log(IAP) for mineral m
double logIAP(int m) const;
//! \brief return log(IAP) for mineral "label mineral"
double logIAP(const std::string& label_mineral) const;
//! \brief return saturation index for mineral m
double saturation_index(int m) const {return logIAP(m) - m_data->logk_mineral(m);}
//! \brief return saturation index for mineral "label mineral"
double saturation_index(const std::string& label_mineral) const;
//! \brief return log(IAP) for mineral m
double logIAP_kinetic(int m) const;
//! \brief return log(IAP) for mineral "label_mineral"
double logIAP_kinetic(const std::string& label_mineral) const;
//! \brief return saturation index for mineral m
double saturation_index_kinetic(int m) const {return logIAP(m) - m_data->logk_mineral_kinetic(m);}
//! \brief return saturation index for mineral m
double saturation_index_kinetic(const std::string& label_mineral) const;
// Gas phase
// =========
//! \brief Return the pressure of gas phase 'g'
double pressure_gas(int g) const {return m_pressure_gas(g);}
//! \brief Return a const reference to the vector of gas pressure
const Eigen::VectorXd& pressure_gas() const {return m_pressure_gas;}
// Thermo variables
// ================
//! \brief Return the temperature
double temperature() const {return get_temperature();}
//! \brief Return the total volume
double volume() const {return get_volume();}
//! \brief Return the volume of the condensed phases
double condensed_phase_volume() const;
//! \brief Return the volume of the gas phase
double gas_phase_volume() const;
// Modification
// ============
double scale_condensed(double scale);
// Database access
// ===============
std::shared_ptr<database::DataContainer> get_database() {return m_data;}
private:
bool m_has_been_properly_initialised;
Eigen::VectorXd m_main_variables;
Eigen::VectorXd m_concentration_secondary;
Eigen::VectorXd m_pressure_gas;
double m_ionic_strength;
Eigen::VectorXd m_loggamma;
std::shared_ptr<database::DataContainer> m_data;
};
} // end namespace specmicp
#undef POW10
#endif // SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP
c4science · Help