Page MenuHomec4science
Files F62535346

adimensional_system_solution_extractor.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 13, 20:50

adimensional_system_solution_extractor.cpp
View Options

#include "adimensional_system_solution_extractor.hpp"
#include "physics/laws.hpp"
#include "database/database.hpp"
namespace specmicp {
scalar_t AdimensionalSystemSolutionExtractor::density_water() const
{
return laws::density_water(units::celsius(25.0), length_unit(), mass_unit());
}
scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_water() const
{
return density_water()*total_saturation_water();
}
scalar_t AdimensionalSystemSolutionExtractor::pH() const
{
// find species responsible for pH
specmicp::database::Database db_handler(m_data);
index_t id = db_handler.component_label_to_id("HO[-]");
if (id != no_species)
{
return 14+log_activity_component(id);
}
else
{
id = db_handler.component_label_to_id("H[+]");
if (id != no_species)
return -log_activity_component(id);
throw std::runtime_error("No component corresponding to the dissociation of water !");
}
}
scalar_t AdimensionalSystemSolutionExtractor::total_saturation_minerals() const
{
return m_solution.main_variables.segment(m_data->nb_component, m_data->nb_mineral).sum();
}
scalar_t AdimensionalSystemSolutionExtractor::mole_concentration_mineral(index_t mineral) const
{
return total_saturation_mineral(mineral)/m_data->molar_volume_mineral(mineral, length_unit());
}
scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_mineral(index_t mineral) const
{
return mole_concentration_mineral(mineral)*m_data->molar_mass_mineral(mineral, mass_unit());
}
//! \brief Return the total aqueous concentration
scalar_t AdimensionalSystemSolutionExtractor::total_aqueous_concentration(index_t component) const
{
scalar_t conc = molality_component(component);
for (index_t aqueous: m_data->range_aqueous())
{
if (m_data->nu_aqueous(aqueous, component) != 0.0)
{
conc += m_data->nu_aqueous(aqueous, component)*molality_aqueous(aqueous);
}
}
return conc;
}
//! \brief Return the total solid concentration
scalar_t AdimensionalSystemSolutionExtractor::total_solid_concentration(index_t component) const
{
scalar_t conc= 0;
for (index_t mineral: m_data->range_mineral())
{
if (m_data->nu_mineral(mineral, component) != 0.0)
{
conc += m_data->nu_mineral(mineral, component)
* total_saturation_mineral(mineral)
/ m_data->molar_volume_mineral(mineral, length_unit());
}
}
return conc;
}
//! \brief Return the total immobile concentration
scalar_t AdimensionalSystemSolutionExtractor::total_immobile_concentration(index_t component) const
{
scalar_t conc= 0;
for (index_t mineral: m_data->range_mineral())
{
if (m_data->nu_mineral(mineral, component) != 0.0)
{
conc += m_data->nu_mineral(mineral, component)
* total_saturation_mineral(mineral)
/ m_data->molar_volume_mineral(mineral, length_unit());
}
}
for (index_t sorbed: m_data->range_sorbed())
{
if (m_data->nu_sorbed(sorbed, component) != 0.0)
{
conc += saturation_water()*density_water()
* m_data->nu_sorbed(sorbed, component)
* sorbed_species_molalities(sorbed);
}
}
return conc;
}
//! Return the saturation index for 'mineral'
scalar_t AdimensionalSystemSolutionExtractor::saturation_index(index_t mineral) const
{
scalar_t saturation_index = - m_data->logk_mineral(mineral);
for (index_t component: m_data->range_aqueous_component())
{
saturation_index += m_data->nu_mineral(mineral, component) *
log_activity_component(component);
}
return saturation_index;
}
//! Return the saturation index for 'mineral_kinetic'
scalar_t AdimensionalSystemSolutionExtractor::saturation_index_kinetic(index_t mineral_kinetic) const
{
scalar_t saturation_index = - m_data->logk_mineral_kinetic(mineral_kinetic);
for (index_t component: m_data->range_aqueous_component())
{
saturation_index += m_data->nu_mineral_kinetic(mineral_kinetic, component) *
log_activity_component(component);
}
return saturation_index;
}
// ########################### //
// //
// Modificator //
// //
// ########################### //
void AdimensionalSystemSolutionModificator::scale_total_concentration(
index_t component,
scalar_t new_value)
{
const scalar_t old_value = total_solid_concentration(component);
const scalar_t factor = new_value/old_value;
m_nonconst_solution.main_variables.segment(dof_mineral(0), m_data->nb_mineral) *= factor;
}
void AdimensionalSystemSolutionModificator::remove_solids()
{
m_nonconst_solution.main_variables.segment(dof_mineral(0), m_data->nb_mineral).setZero();
}
Vector AdimensionalSystemSolutionModificator::set_minerals_kinetics(std::vector<index_t>& list_species)
{
index_t nb_kinetics = list_species.size();
index_t nb_new_mineral = m_data->nb_mineral - nb_kinetics;
std::vector<index_t> minerals_to_keep;
minerals_to_keep.reserve(nb_new_mineral);
std::vector<index_t> new_kinetics_index(nb_kinetics, no_species);
Vector saturation_kinetics(nb_kinetics);
index_t new_ind_eq = m_data->nb_component;
index_t new_ind_kin = 0;
index_t tot_ind_kin = m_data->nb_mineral_kinetic;
// ###TODO optimize
for (index_t mineral: m_data->range_mineral())
{
auto is_kin = std::find(list_species.begin(), list_species.end(), mineral);
// If mineral is still at equilibrium
if (is_kin == list_species.end())
{
minerals_to_keep.push_back(mineral);
m_nonconst_solution.main_variables(new_ind_eq) = total_saturation_mineral(mineral);
++new_ind_eq;
}
// If mineral is governed by kinetics
else
{
saturation_kinetics(new_ind_kin) = total_saturation_mineral(mineral);
++new_ind_kin;
// save the new index (index in the kinetics vector)
// The order is conserved !
new_kinetics_index[is_kin - list_species.begin()] = tot_ind_kin;
++tot_ind_kin;
}
}
specmicp_assert(new_ind_eq == m_data->nb_component+nb_new_mineral);
// change the database
database::Database dbhandler(m_data);
dbhandler.minerals_keep_only(minerals_to_keep);
// update the list of species
list_species.swap(new_kinetics_index);
// resize
m_nonconst_solution.main_variables.conservativeResize(m_data->nb_component+nb_new_mineral);
return saturation_kinetics;
}
} // end namespace specmicp

Event Timeline

c4science · Help