reactant_box.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Oct 18, 15:00

reactant_box.cpp
View Options

	/* =============================================================================

	Copyright (c) 2014 - 2016
	F. 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:

	1. Redistributions of source code must retain the above copyright notice,
	this list of conditions and the following disclaimer.

	2. 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.

	3. Neither the name of the copyright holder 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 HOLDER 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 "reactant_box.hpp"

	#include "specmicp/adimensional/adimensional_system_structs.hpp"

	#include "specmicp_database/data_container.hpp"
	#include "specmicp_database/database.hpp"
	#include "specmicp_database/unknown_class.hpp"

	#include "specmicp_common/physics/laws.hpp"

	#include <vector>
	#include <unordered_map>
	#include <stdexcept>

	#include <iostream>

	namespace specmicp {


	struct SPECMICP_DLL_LOCAL fug_container {
	std::string gas;
	std::string component;
	scalar_t fugacity;

	fug_container(
	const std::string& agas,
	const std::string& acomponent,
	scalar_t t_fugacity
	):
	gas(agas), component(acomponent), fugacity(t_fugacity)
	{}
	};

	struct ReactantBox::ReactantBoxImpl {
	scalar_t mass_solution {-1};
	std::vector<std::string> keep_components {};
	std::unordered_map<std::string, scalar_t> aqueous_species {};
	std::unordered_map<std::string, scalar_t> solid_phase {};

	// constraints
	std::string charge_keeper {""};
	scalar_t inert_vol_frac {0.0};
	WaterEquationType water_eq_type {WaterEquationType::MassConservation};
	scalar_t water_param {-1.0};
	water_partial_pressure_f h2o_partial_pressure_model;
	std::vector<fug_container> fix_fugacity;
	std::vector<std::pair<std::string, scalar_t>> fix_activity;
	std::vector<std::pair<std::string, scalar_t>> fix_molality;
	std::vector<std::pair<std::string, scalar_t>> mineral_upper_bound;

	// database lookup
	database::SpeciesTypeFinder db_search;

	ReactantBoxImpl(RawDatabasePtr thedb):
	db_search(thedb)
	{}

	std::vector<index_t> get_components_to_keep() const;

	void dissolve_aqueous_species(
	Vector& total_concentration
	) const;
	void dissolve_solid_phases(
	Vector& total_concentration
	) const;

	void raise_unknown_unit(const std::string& error_msg) const {
	throw std::invalid_argument(error_msg);
	}
	void assert_mass_solution() const {
	if (mass_solution == -1) {
	throw std::logic_error("Mass of solution is not set !");
	}
	}

	scalar_t factor_total_concentration(const units::UnitsSet& units_set) const;

	scalar_t parse_unit_aqueous_species(
	std::string name,
	const units::AmountUnit& unit
	) const;
	scalar_t parse_unit_solid_phase(
	std::string name,
	const units::AmountUnit& unit
	) const;
	};

	ReactantBox::ReactantBox(
	std::shared_ptr<database::DataContainer> raw_data,
	const units::UnitsSet& units_set
	):
	units::UnitBaseClass(units_set),
	database::DatabaseHolder(raw_data),
	m_impl(utils::make_pimpl<ReactantBoxImpl>(raw_data))
	{

	}


	ReactantBox::ReactantBox(ReactantBox&& other):
	units::UnitBaseClass(other.get_units()),
	database::DatabaseHolder(other.m_data),
	m_impl(std::move(other.m_impl))
	{

	}


	ReactantBox& ReactantBox::operator=(ReactantBox&& other)
	{
	m_impl = std::move(other.m_impl);
	return *this;
	}

	ReactantBox::~ReactantBox() = default;

	void ReactantBox::set_solution(scalar_t value, std::string str_unit)
	{
	set_solution(value, units::parse_amount_unit(str_unit));
	}
	void ReactantBox::set_solution(
	scalar_t value,
	const units::AmountUnit& unit
	)
	{
	scalar_t factor {1.0};
	if (unit.type == units::AmountUnitType::Mass) {
	factor = unit.factor_si;
	} else if (unit.type == units::AmountUnitType::Volume) {
	factor = unit.factor_si*laws::density_water(units::SI_units);
	} else if (unit.type == units::AmountUnitType::VolumeFraction) {
	factor = unit.factor_si*laws::density_water(units::SI_units);
	} else if (unit.type == units::AmountUnitType::MassConcentration) {
	factor = unit.factor_si;
	} else {
	m_impl->raise_unknown_unit(
	"Only mass. volume, volume fraction,"
	" or mass concentration units "
	"are accepted for the amount of solution"
	);
	}
	m_impl->mass_solution = factor*value;

	}

	void ReactantBox::add_aqueous_species(
	std::string name,
	scalar_t value,
	std::string str_unit
	)
	{
	add_aqueous_species(name, value, units::parse_amount_unit(str_unit));
	}

	void ReactantBox::add_aqueous_species(
	std::string name,
	scalar_t value,
	const units::AmountUnit& unit
	)
	{

	const scalar_t factor = m_impl->parse_unit_aqueous_species(name, unit);
	const scalar_t molality = factor*value;
	m_impl->aqueous_species[name] += molality;
	}

	void ReactantBox::set_aqueous_species(
	std::string name,
	scalar_t value,
	std::string str_unit
	)
	{
	set_aqueous_species(name, value, units::parse_amount_unit(str_unit));
	}

	void ReactantBox::set_aqueous_species(
	std::string name,
	scalar_t value,
	const units::AmountUnit& unit
	)
	{

	const scalar_t factor = m_impl->parse_unit_aqueous_species(name, unit);
	const scalar_t molality = factor*value;
	m_impl->aqueous_species[name] = molality;
	}

	void ReactantBox::add_solid_phase(
	std::string name,
	scalar_t value,
	std::string str_unit
	)
	{
	return add_solid_phase(name, value, units::parse_amount_unit(str_unit));
	}

	void ReactantBox::add_solid_phase(
	std::string name,
	scalar_t value,
	const units::AmountUnit& unit
	)
	{

	const scalar_t factor = m_impl->parse_unit_solid_phase(name, unit);
	m_impl->solid_phase[name] += factor*value;
	}

	void ReactantBox::set_solid_phase(
	std::string name,
	scalar_t value,
	std::string str_unit
	)
	{
	return set_solid_phase(name, value, units::parse_amount_unit(str_unit));
	}

	void ReactantBox::set_solid_phase(
	std::string name,
	scalar_t value,
	const units::AmountUnit& unit
	)
	{

	const scalar_t factor = m_impl->parse_unit_solid_phase(name, unit);
	m_impl->solid_phase[name] = factor*value;
	}

	void ReactantBox::add_component(const std::string& name)
	{
	m_impl->keep_components.push_back(name);
	}

	Vector ReactantBox::get_total_concentration(bool modify_db) const
	{
	std::shared_ptr<database::DataContainer> raw_data = get_database();

	if (modify_db) {
	// remove extra components
	const auto id_comp_to_keep = m_impl->get_components_to_keep();
	database::Database db_manager(raw_data);
	db_manager.keep_only_components(id_comp_to_keep);
	}


	Vector total_concentration = Vector::Zero(raw_data->nb_component());
	total_concentration[database::water_id] =
	m_impl->mass_solution/raw_data->molar_mass_basis(database::water_id, units::SI_units);
	m_impl->dissolve_aqueous_species(total_concentration);
	m_impl->dissolve_solid_phases(total_concentration);

	total_concentration *= m_impl->factor_total_concentration(get_units());

	return total_concentration;
	}

	void ReactantBox::ReactantBoxImpl::dissolve_aqueous_species(
	Vector& total_concentration
	) const
	{
	auto raw_data = db_search.get_database();
	for (auto& it: aqueous_species) {
	database::GenericAqueousSpecies aq = db_search.find_aqueous_species(it.first);
	const scalar_t conc = mass_solution*it.second;
	if (aq.type == database::AqueousSpeciesClass::Component)
	{
	total_concentration[aq.id] += conc;
	}
	else if (aq.type == database::AqueousSpeciesClass::Aqueous)
	{
	for (auto idc: raw_data->range_component())
	{
	const auto nu = raw_data->nu_aqueous(aq.id, idc);
	if (nu != 0)
	{
	total_concentration[idc] += nu*conc;
	}
	}
	}
	else if (aq.type == database::AqueousSpeciesClass::Compound)
	{
	for (auto idc: raw_data->range_component())
	{
	const auto nu = raw_data->nu_compound(aq.id, idc);
	if (nu != 0)
	{
	total_concentration[idc] += nu*conc;
	}
	}
	}
	else
	{
	throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species.");
	}
	}
	}


	void ReactantBox::ReactantBoxImpl::dissolve_solid_phases(
	Vector& total_concentration
	) const
	{
	auto raw_data = db_search.get_database();
	for (auto& it: solid_phase) {
	database::GenericSolidPhase sp = db_search.find_solid_phase(it.first);
	if (sp.type == database::SolidPhaseClass::EquilibriumMineral)
	{
	for (auto idc: raw_data->range_component()) {
	const auto nu = raw_data->nu_mineral(sp.id, idc);
	if (nu != 0) {
	total_concentration[idc] += nu*it.second;
	}
	}
	}
	else if (sp.type == database::SolidPhaseClass::MineralKinetics)
	{
	for (auto idc: raw_data->range_component()) {
	const auto nu = raw_data->nu_mineral_kinetic(sp.id, idc);
	if (nu != 0) {
	total_concentration[idc] += nu*it.second;
	}
	}
	}
	else
	{
	throw std::invalid_argument("'"+it.first+"' is not a valid solid phase.");
	}
	}
	}

	std::vector<index_t> ReactantBox::ReactantBoxImpl::get_components_to_keep(
	) const
	{
	auto raw_data = db_search.get_database();

	std::vector<bool> to_keep(raw_data->nb_component(), false);
	to_keep[database::water_id] = true;
	to_keep[database::electron_id] = true;

	// aqueous species
	// ---------------
	for (auto& it: aqueous_species) {
	database::GenericAqueousSpecies aq = db_search.find_aqueous_species(it.first);
	if (aq.type == database::AqueousSpeciesClass::Component)
	{
	to_keep[aq.id] = true;
	}
	else if (aq.type == database::AqueousSpeciesClass::Aqueous)
	{
	for (auto idc: raw_data->range_aqueous_component())
	{
	if (raw_data->nu_aqueous(aq.id, idc) != 0)
	{
	to_keep[idc] = true;
	}
	}
	}
	else if (aq.type == database::AqueousSpeciesClass::Compound)
	{
	for (auto idc: raw_data->range_aqueous_component())
	{
	if (raw_data->nu_compound(aq.id, idc) != 0)
	{
	to_keep[idc] = true;
	}
	}
	}
	else
	{
	throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species.");
	}
	}
	// solid phases
	// ------------
	for (auto& it: solid_phase) {
	database::GenericSolidPhase sp = db_search.find_solid_phase(it.first);
	if (sp.type == database::SolidPhaseClass::EquilibriumMineral)
	{
	for (auto idc: raw_data->range_aqueous_component())
	{
	if (raw_data->nu_mineral(sp.id, idc) != 0)
	{
	to_keep[idc] = true;
	}
	}
	}
	else if (sp.type == database::SolidPhaseClass::MineralKinetics)
	{
	for (auto idc: raw_data->range_aqueous_component())
	{
	if (raw_data->nu_mineral_kinetic(sp.id, idc) != 0)
	{
	to_keep[idc] = true;
	}
	}
	}
	else
	{
	throw std::invalid_argument("'"+it.first+"' is not a valid solid phase.");
	}
	}
	// extra components
	// -----------------
	for (auto& comp_label: keep_components) {
	const index_t idc = raw_data->get_id_component(comp_label);
	specmicp_assert(idc != no_species);
	to_keep[idc] = true;
	}

	// get ids of components to keep
	std::vector<index_t> id_comp_to_keep;
	id_comp_to_keep.reserve(raw_data->nb_component());

	for (auto id: raw_data->range_component()) {
	if (to_keep[id]) {
	id_comp_to_keep.push_back(id);
	}
	}
	return id_comp_to_keep;

	}

	scalar_t ReactantBox::ReactantBoxImpl::factor_total_concentration(
	const units::UnitsSet& units_set
	) const
	{
	scalar_t mol_fac = units::scaling_factor(units_set.quantity);
	scalar_t length_fac = units::scaling_factor(units_set.length);

	return std::pow(length_fac, 3)/mol_fac;
	}


	void ReactantBox::set_charge_keeper(std::string charge_keeper)
	{
	auto id = m_data->get_id_component(charge_keeper);
	if (id == no_species) // check that the component exist
	{
	// if not check if it is an element and find the corresponding component
	id = m_data->get_id_component_from_element(charge_keeper);
	if (id == no_species)
	{
	throw std::invalid_argument(
	"ReactantBox::set_charge_keeper : "
	"the argument (" + charge_keeper
	+ ") is not a component nor an element");
	}
	charge_keeper = m_data->get_label_component(id);
	}

	m_impl->charge_keeper = charge_keeper;
	add_component(charge_keeper); // so it is always in the system
	}
	void ReactantBox::set_inert_volume_fraction(scalar_t inert_volume_fraction)
	{
	m_impl->inert_vol_frac = inert_volume_fraction;
	}
	void ReactantBox::set_saturated_system()
	{
	m_impl->water_eq_type = WaterEquationType::SaturatedSystem;
	}
	void ReactantBox::set_fixed_saturation(scalar_t saturation)
	{
	if (saturation <= 0.0 or saturation >= 1.0) {
	throw std::invalid_argument("The saturation must be between 0 and 1"
	"(Value : " + std::to_string(saturation) + ")");
	}
	m_impl->water_eq_type = WaterEquationType::FixedSaturation;
	m_impl->water_param = saturation;
	}
	void ReactantBox::disable_conservation_water()
	{
	m_impl->water_eq_type = WaterEquationType::NoEquation;
	}

	void ReactantBox::set_water_partial_pressure_model(
	water_partial_pressure_f h2o_pressure_model
	)
	{
	m_impl->h2o_partial_pressure_model = h2o_pressure_model;
	}
	void ReactantBox::add_fixed_fugacity_gas(
	std::string gas,
	std::string component,
	scalar_t fugacity
	)
	{
	specmicp_assert(fugacity > 0);
	m_impl->fix_fugacity.emplace_back(gas, component, fugacity);
	add_component(component);
	}
	void ReactantBox::add_fixed_activity_component(
	std::string component,
	scalar_t activity
	)
	{
	specmicp_assert(activity > 0);
	m_impl->fix_activity.emplace_back(component, activity);
	add_component(component);
	}
	void ReactantBox::add_fixed_molality_component(
	std::string component,
	scalar_t molality
	)
	{
	specmicp_assert(molality > 0);
	m_impl->fix_molality.emplace_back(component, molality);
	add_component(component);
	}

	void ReactantBox::set_mineral_upper_bound(
	std::string mineral,
	scalar_t max_volume_fraction
	)
	{
	if (max_volume_fraction < 0 or (not std::isfinite(max_volume_fraction))) {
	throw std::invalid_argument("The saturation must a finite positive number"
	"(Value : " + std::to_string(max_volume_fraction) + ")");
	}
	m_impl->mineral_upper_bound.emplace_back(mineral, max_volume_fraction);
	}


	AdimensionalSystemConstraints ReactantBox::get_constraints(bool modify_db) const
	{
	database::DataContainer* raw_data = get_database().get();

	AdimensionalSystemConstraints constraints;
	constraints.total_concentrations = get_total_concentration(modify_db);

	// Water
	switch (m_impl->water_eq_type) {
	case WaterEquationType::MassConservation:
	constraints.enable_conservation_water();
	break;
	case WaterEquationType::FixedSaturation:
	specmicp_assert(m_impl->water_param > 0.0 and m_impl->water_param < 1.0);
	constraints.set_fixed_saturation(m_impl->water_param);
	break;
	case WaterEquationType::SaturatedSystem:
	constraints.set_saturated_system();
	break;
	case WaterEquationType::NoEquation:
	constraints.disable_conservation_water();
	}
	if (m_impl->h2o_partial_pressure_model != nullptr)
	{
	constraints.set_water_partial_pressure_model(m_impl->h2o_partial_pressure_model);
	}

	// volume fraction
	constraints.set_inert_volume_fraction(
	m_impl->inert_vol_frac
	);


	// Components

	if (m_impl->charge_keeper != "")
	{
	constraints.set_charge_keeper(
	raw_data->get_id_component(m_impl->charge_keeper)
	);
	}

	for (auto& fug: m_impl->fix_fugacity)
	{
	constraints.add_fixed_fugacity_gas(
	raw_data->get_id_gas(fug.gas),
	raw_data->get_id_component(fug.component),
	std::log10(fug.fugacity)
	);
	}
	for (auto& act: m_impl->fix_activity)
	{
	constraints.add_fixed_activity_component(
	raw_data->get_id_component(act.first),
	std::log10(act.second)
	);
	}
	for (auto& mol: m_impl->fix_molality)
	{
	constraints.add_fixed_molality_component(
	raw_data->get_id_component(mol.first),
	std::log10(mol.second)
	);
	}

	for (auto& mup: m_impl->mineral_upper_bound)
	{
	constraints.set_mineral_upper_bound(
	raw_data->get_id_mineral(mup.first),
	mup.second
	);
	}

	return constraints;
	}



	scalar_t ReactantBox::ReactantBoxImpl::parse_unit_aqueous_species(
	std::string name,
	const units::AmountUnit& unit
	) const
	{
	scalar_t factor = 1.0;
	if (unit.type == units::AmountUnitType::Molality)
	{
	factor = unit.factor_si;
	}
	else if (unit.type == units::AmountUnitType::MoleConcentration)
	{
	factor = unit.factor_si/laws::density_water(units::SI_units);
	}
	else if (unit.type == units::AmountUnitType::MassConcentration)
	{
	const scalar_t molar_mass = db_search.molar_mass_aqueous(name);
	factor = unit.factor_si / molar_mass
	/ laws::density_water(units::SI_units);
	}
	else if (unit.type == units::AmountUnitType::Mass)
	{
	assert_mass_solution();
	const scalar_t molar_mass = db_search.molar_mass_aqueous(name);
	factor = unit.factor_si / molar_mass / mass_solution;
	}
	else if (unit.type == units::AmountUnitType::NumberOfMoles)
	{
	assert_mass_solution();
	factor = unit.factor_si / mass_solution;
	}
	else
	{
	raise_unknown_unit("Unit for aqueous species " + name + "is invalid."
	" Only accepted units for aqueous species are"
	" molality, mole and mass concettration,"
	" mass and number of moles.");
	}
	return factor;
	}

	scalar_t ReactantBox::ReactantBoxImpl::parse_unit_solid_phase(
	std::string name,
	const units::AmountUnit& unit
	) const
	{
	scalar_t factor = 1.0;

	// convert to mole concentration is SI
	if (unit.type == units::AmountUnitType::MoleConcentration or
	unit.type == units::AmountUnitType::NumberOfMoles)
	{
	factor = unit.factor_si;
	}
	else if (unit.type == units::AmountUnitType::MassConcentration or
	unit.type == units::AmountUnitType::Mass
	)
	{
	const scalar_t molar_mass = db_search.molar_mass_solid_phase(name);
	factor = unit.factor_si / molar_mass;
	}
	else if (unit.type == units::AmountUnitType::Volume)
	{
	const scalar_t molar_volume = db_search.molar_volume_solid_phase(name);
	factor = unit.factor_si / molar_volume;
	}
	else if (unit.type == units::AmountUnitType::VolumeFraction)
	{
	const scalar_t molar_volume = db_search.molar_volume_solid_phase(name);
	factor = unit.factor_si / molar_volume;

	}
	else
	{
	raise_unknown_unit("Unit for solid phase " + name + " is invalid.");
	}
	return factor;
	}

	} // end namespace specmicp

reactant_box.cppNo OneTemporaryActions

File Metadata

reactant_box.cppView Options

Event Timeline

reactant_box.cpp
No OneTemporary
Actions

reactant_box.cpp
View Options