Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/equilibrium_data.hpp
b9a808db21fbreactivetransport

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 "common.hpp"
#include "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,
scalar_t 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 Vector& main_vars,
const Vector& secondary,
const Vector& loggamma,
scalar_t 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 Vector& main_vars,
const Vector& secondary,
const Vector& loggamma,
scalar_t ionic_strength,
const Vector& pressures,
scalar_t 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
scalar_t mass_water() const {return m_main_variables(0);}
//! \brief Return the molality of component i
scalar_t molality_component(index_t i) const {return POW10(m_main_variables(i));}
//! \brief Return log10 of the molality of component i
scalar_t log_molality_component(index_t i) const {return m_main_variables(i);}
//! \brief Return the number of moles of the mineral at equilibrium 'm'
scalar_t moles_mineral(index_t 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)
scalar_t mass_mineral(index_t m) const;
//! \brief Return the volume of a mineral (V)
scalar_t volume_mineral(index_t m) const {
return m_data->molar_volume_mineral(m)*moles_mineral(m);
}
//! \brief Return the volume of all minerals
scalar_t moles_minerals() const;
//! \brief Return the volume of all minerals
scalar_t volume_minerals() const;
//! \brief Return a const reference to the vector of main variables
const Vector& main_variables() const {return m_main_variables;}
//! \brief Return the total concentrations
void total_concentrations(Vector& total_conc) const;
// molality
// ========
//! \brief return the pH
scalar_t pH() const;
//! \brief Return the molality of secondary species j
scalar_t 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 Vector& molality_secondary() const {return m_concentration_secondary;}
//! \brief Return total concentration of aqueous species for component i
scalar_t total_aqueous_concentration_component(index_t i) const;
//! \brief Compute total aqeuous concentration for each component
void total_aqueous_concentrations(Vector& totconc) const;
// Activity
// ========
scalar_t ionic_strength() const {return m_ionic_strength;}
//! \brief Return a const reference to the activity coefficient vector
const Vector& activity_coefficient() const {return m_loggamma;}
//! \brief Return activity coefficient of component i
scalar_t activity_coefficient_component(int i) const {
assert(i >= 1);
return POW10(m_loggamma(i));
}
//! \brief Return logarithm of activity coefficient of component i
scalar_t 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
scalar_t 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
scalar_t 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)
scalar_t activity_component(int i) const {
assert(i >= 1);
return POW10(m_loggamma(i) + m_main_variables(i));
}
//! \brief Return activity of secondary species j
scalar_t 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
// ========
void remove_minerals() {
m_main_variables.segment(m_data->nb_component, m_data->nb_mineral).setZero();}
//! \brief return log(IAP) for mineral m
scalar_t logIAP(index_t m) const;
//! \brief return log(IAP) for mineral "label mineral"
scalar_t logIAP(const std::string& label_mineral) const;
//! \brief return saturation index for mineral m
scalar_t saturation_index(index_t m) const {return logIAP(m) - m_data->logk_mineral(m);}
//! \brief return saturation index for mineral "label mineral"
scalar_t saturation_index(const std::string& label_mineral) const;
//! \brief return log(IAP) for mineral m
scalar_t logIAP_kinetic(index_t m) const;
//! \brief return log(IAP) for mineral "label_mineral"
scalar_t logIAP_kinetic(const std::string& label_mineral) const;
//! \brief return saturation index for mineral m
scalar_t saturation_index_kinetic(int m) const {return logIAP(m) - m_data->logk_mineral_kinetic(m);}
//! \brief return saturation index for mineral m
scalar_t saturation_index_kinetic(const std::string& label_mineral) const;
// Gas phase
// =========
//! \brief Return the pressure of gas phase 'g'
scalar_t pressure_gas(int g) const {return m_pressure_gas(g);}
//! \brief Return a const reference to the vector of gas pressure
const Vector& pressure_gas() const {return m_pressure_gas;}
// Thermo variables
// ================
//! \brief Return the temperature
scalar_t temperature() const {return get_temperature();}
//! \brief Return the total volume
scalar_t volume() const {return get_volume();}
//! \brief Return the volume of the condensed phases
scalar_t condensed_phase_volume() const;
//! \brief Return the volume of the gas phase
scalar_t gas_phase_volume() const;
// Modification
// ============
void scale_condensed(scalar_t scale);
// Database access
// ===============
std::shared_ptr<database::DataContainer> get_database() const {return m_data;}
private:
bool m_has_been_properly_initialised;
Vector m_main_variables;
Vector m_concentration_secondary;
Vector m_pressure_gas;
scalar_t m_ionic_strength;
Vector m_loggamma;
database::RawDatabasePtr m_data;
};
} // end namespace specmicp
#undef POW10
#endif // SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP
c4science · Help