diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp
index afa2d15..2796376 100644
--- a/src/specmicp/adimensional/adimensional_system.cpp
+++ b/src/specmicp/adimensional/adimensional_system.cpp
@@ -1,1878 +1,1881 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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 "adimensional_system.hpp"
#include "adimensional_system_solution.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_common/physics/constants.hpp"
#include "specmicp_common/physics/laws.hpp"
#include <cmath>
#include <random>
//#include <iostream>
// uncomment to activate the finite difference jacobian
// deprecated: use the cmake option instead
// #define SPECMICP_DEBUG_EQUATION_FD_JACOBIAN
namespace specmicp {
static const scalar_t log10 = std::log(10.0);
// Constructor
// ===========
// No previous solution
// --------------------
AdimensionalSystem::AdimensionalSystem(
RawDatabasePtr& ptrdata,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemOptions& options,
const units::UnitsSet& units_set
):
AdimemsionalSystemNumbering(ptrdata),
OptionsHandler<AdimensionalSystemOptions>(options),
units::UnitBaseClass(units_set),
m_inert_volume_fraction(constraints.inert_volume_fraction),
m_second(ptrdata),
m_equations(total_dofs(), ptrdata)
{
specmicp_assert(ptrdata->is_valid());
m_fixed_values.setZero(ptrdata->nb_component()+ptrdata->nb_ssites()+1);
number_eq(constraints);
}
// Previous solution
// -----------------
AdimensionalSystem::AdimensionalSystem(
RawDatabasePtr& ptrdata,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolution& previous_solution,
const AdimensionalSystemOptions& options,
const units::UnitsSet& units_set
):
AdimemsionalSystemNumbering(ptrdata),
OptionsHandler<AdimensionalSystemOptions>(options),
units::UnitBaseClass(units_set),
m_inert_volume_fraction(constraints.inert_volume_fraction),
m_second(previous_solution),
m_equations(total_dofs(), ptrdata)
{
specmicp_assert(ptrdata->is_valid());
m_fixed_values.setZero(ptrdata->nb_component()+ptrdata->nb_ssites()+1);
number_eq(constraints);
}
// Secondary variables constructor
// ===============================
// No previous solution
// --------------------
AdimensionalSystem::SecondaryVariables::SecondaryVariables(
const RawDatabasePtr& data
):
secondary_molalities(Vector::Zero(data->nb_aqueous())),
loggamma(Vector::Zero(data->nb_component()+data->nb_aqueous())),
gas_fugacity(Vector::Zero(data->nb_gas())),
gas_concentration(Vector::Zero(data->nb_gas())),
sorbed_concentrations(Vector::Zero(data->nb_sorbed()))
{}
// Previous solution
// -----------------
AdimensionalSystem::SecondaryVariables::SecondaryVariables(
const AdimensionalSystemSolution& previous_solution
):
secondary_molalities(previous_solution.secondary_molalities),
loggamma(previous_solution.log_gamma),
gas_fugacity(previous_solution.gas_fugacities),
gas_concentration(Vector::Zero(previous_solution.gas_fugacities.rows())),
sorbed_concentrations(previous_solution.sorbed_molalities)
{}
// IdEquations constructor
// =======================
AdimensionalSystem::IdEquations::IdEquations(
index_t nb_dofs,
const RawDatabasePtr& data
):
ideq(nb_dofs, no_equation),
component_equation_type(data->nb_component()+data->nb_ssites()+1, no_equation),
fixed_activity_species(data->nb_component()+1, no_species),
active_aqueous(data->nb_aqueous(), false),
active_gas(data->nb_gas(), false),
active_sorbed(data->nb_sorbed())
{}
// Equation numbering
// ==================
// Note : this function also computes scaling factor that would be used in the computation
// ------
void AdimensionalSystem::number_eq(
const AdimensionalSystemConstraints& constraints
)
{
index_t neq = 0;
// units
compute_scaling_factors();
// Water
// =====
if (constraints.water_equation != WaterEquationType::NoEquation)
{
m_equations.type_equation(dof_water()) = static_cast<index_t>(constraints.water_equation);
if (constraints.water_equation == WaterEquationType::MassConservation)
{
m_fixed_values(dof_water()) = constraints.total_concentrations(dof_water());
if (constraints.water_partial_pressure.use_partial_pressure_model)
{
m_equations.use_water_pressure_model = true;
m_equations.water_pressure_model = constraints.water_partial_pressure.partial_pressure_model;
}
}
else if (constraints.water_equation == WaterEquationType::FixedSaturation)
{
m_fixed_values(dof_water()) = constraints.water_parameter;
}
else
{
specmicp_assert(constraints.water_equation == WaterEquationType::SaturatedSystem);
}
m_equations.add_equation(dof_water(), &neq);
}
// Aqueous components
// ==================
number_eq_aqueous_component(constraints, neq);
// Surface model
// =============
// only activate surface model if at least one concentration is positive
if (constraints.surface_model.model_type == SurfaceEquationType::Equilibrium)
{
specmicp_assert(constraints.surface_model.total_ssite_concentrations.rows() == m_data->nb_ssites());
bool at_least_one = false;
for (index_t q: m_data->range_ssites()) {
const scalar_t tot_conc = constraints.surface_model.total_ssite_concentrations(q);
if (tot_conc > 0.0) {
// add the equation
m_equations.add_equation(dof_surface(q), &neq);
m_equations.type_equation(dof_surface(q)) = static_cast<index_t>(SurfaceEquationType::Equilibrium);
// setup the total concentration
m_fixed_values(dof_surface(q)) = tot_conc;
at_least_one = true;
} else {
m_equations.type_equation(dof_surface(q)) = static_cast<index_t>(SurfaceEquationType::NoEquation);
}
}
// No surface potential equation for this model
if (at_least_one) {
m_equations.type_equation(dof_surface_potential()) = static_cast<index_t>(SurfaceEquationType::Equilibrium);
} else {
m_equations.type_equation(dof_surface_potential()) = static_cast<index_t>(SurfaceEquationType::NoEquation);
}
m_fixed_values(dof_surface_potential()) = -1.0;
}
else if (constraints.surface_model.model_type == SurfaceEquationType::EDL)
{
specmicp_assert(constraints.surface_model.total_ssite_concentrations.rows() == m_data->nb_ssites());
bool at_least_one = false;
for (index_t q: m_data->range_ssites()) {
const scalar_t tot_conc = constraints.surface_model.total_ssite_concentrations(q);
if (tot_conc > 0.0) {
// add the equation
m_equations.add_equation(dof_surface(q), &neq);
m_equations.type_equation(dof_surface(q)) = static_cast<index_t>(SurfaceEquationType::EDL);
// setup the total concentration
m_fixed_values(dof_surface(q)) = tot_conc;
at_least_one = true;
} else {
m_equations.type_equation(dof_surface(q)) = static_cast<index_t>(SurfaceEquationType::NoEquation);
}
}
if (at_least_one) {
m_equations.type_equation(dof_surface_potential()) = static_cast<index_t>(SurfaceEquationType::EDL);
m_equations.add_equation(dof_surface_potential(), &neq);
m_fixed_values(dof_surface_potential()) = constraints.surface_model.sorption_surface;
} else {
m_equations.type_equation(dof_surface_potential()) = static_cast<index_t>(SurfaceEquationType::NoEquation);
}
}
else {
m_equations.type_equation(dof_surface_potential()) = static_cast<index_t>(SurfaceEquationType::NoEquation);
}
// Secondary species
// =================
// Secondary aqueous species
// -------------------------
bool include_half_cell_reaction = (constraints.electron_constraint.equation_type != ElectronEquationType::NoEquation);
bool solve_electron_equation {false};
for (auto j: m_data->range_aqueous())
{
bool can_exist { true };
if ( include_half_cell_reaction or not m_data->is_half_cell_reaction(j))
for (const auto& k: m_equations.nonactive_component)
{
if (m_data->nu_aqueous(j, k) != 0.0)
{
can_exist = false;
break;
}
}
else
{
can_exist = false;
}
m_equations.set_aqueous_active(j, can_exist);
if (can_exist and m_data->is_half_cell_reaction(j))
solve_electron_equation = true;
//std::cout << m_data->labels_aqueous[j] << "can exist ? " << can_exist << " - logk : " << m_data->logk_aqueous(j) << std::endl;
}
// Gas
// ---
for (index_t k: m_data->range_gas())
{
bool can_exist = true;
for (const index_t& n: m_equations.nonactive_component)
{
if (m_data->nu_gas(k, n) != 0.0)
{
can_exist = false;
break;
}
}
m_equations.set_gas_active(k, can_exist);
}
// Sorbed species
// --------------
- if (surface_model() == SurfaceEquationType::NoEquation) {
+ if (surface_model() == SurfaceEquationType::NoEquation or constraints.immobile_disabled) {
for (index_t s: m_data->range_sorbed()) {
m_equations.set_sorbed_active(s, false);
}
}
else {
// non active sorbed : no surface site or non active component
for (index_t s: m_data->range_sorbed()) {
if (surface_total_concentration(m_data->sorbed.surface_site(s)) == 0) {
m_equations.set_sorbed_active(s, false);
} else {
bool can_exist = true;
for (const index_t& k: m_equations.nonactive_component)
{
if (m_data->nu_sorbed(s, k) != 0.0)
{
can_exist = false;
break;
}
}
m_equations.set_sorbed_active(s, can_exist);
}
}
}
// Electron equation
// -----------------
if (solve_electron_equation)
{
m_equations.add_equation(dof_electron(), &neq);
m_equations.type_equation(dof_electron()) =
static_cast<index_t>(constraints.electron_constraint.equation_type);
if (constraints.electron_constraint.equation_type == ElectronEquationType::Equilibrium)
{
m_fixed_values(dof_electron()) = 0.0;
}
else if (constraints.electron_constraint.equation_type == ElectronEquationType::FixedpE)
{
m_fixed_values(dof_electron()) = constraints.electron_constraint.fixed_value;
m_equations.related_species(dof_electron()) = constraints.electron_constraint.species;
//assert(m_fixed_activity_species[dof_electron()] >= 0
// and m_fixed_activity_species[dof_electron()] < m_data->nb_aqueous());
//assert(m_data->is_half_cell_reaction(m_fixed_activity_species[dof_electron()]));
}
// scaling
if (get_options().scaling_electron == 0.0)
{
for (index_t component : m_data->range_aqueous_component())
{
if (aqueous_component_equation_type(component)
== AqueousComponentEquationType::MassConservation)
{
get_options().scaling_electron = total_concentration_bc(component);
break;
}
}
}
}
// above equations are 'free' (i.e. non constrained)
m_equations.nb_free_variables = neq;
// following equations are complementarity conditions
// Minerals
// ========
- for (index_t m: m_data->range_mineral())
+ if (not constraints.immobile_disabled)
{
- bool can_precipitate = true; // true if all the components of the mineral
- // are present in the system
- // just check that the molar volume exist
- auto molar_volume = m_data->molar_volume_mineral(m);
- // Remove minerals that cannot precipitate
- for (index_t& k: m_equations.nonactive_component)
+ for (index_t m: m_data->range_mineral())
{
- if (m_data->nu_mineral(m, k) != 0.0 and molar_volume > 0.0)
+ bool can_precipitate = true; // true if all the components of the mineral
+ // are present in the system
+ // just check that the molar volume exist
+ auto molar_volume = m_data->molar_volume_mineral(m);
+ // Remove minerals that cannot precipitate
+ for (index_t& k: m_equations.nonactive_component)
{
- can_precipitate = false;
- break; // this is not a mineral that can precipitate
+ if (m_data->nu_mineral(m, k) != 0.0 and molar_volume > 0.0)
+ {
+ can_precipitate = false;
+ break; // this is not a mineral that can precipitate
+ }
}
- }
- // check if no precipitation is set
- // if upper bound is 0, deactivate the mineral
- if (not constraints.mineral_constraints.empty())
- {
- for (const MineralConstraint& v: constraints.mineral_constraints)
+ // check if no precipitation is set
+ // if upper bound is 0, deactivate the mineral
+ if (not constraints.mineral_constraints.empty())
{
- if (v.id_mineral == m)
+ for (const MineralConstraint& v: constraints.mineral_constraints)
{
- if (v.param == 0) {
- can_precipitate = false;
+ if (v.id_mineral == m)
+ {
+ if (v.param == 0) {
+ can_precipitate = false;
+ }
}
}
}
- }
- if (can_precipitate)
- {
- m_equations.add_equation(dof_mineral(m), &neq);
+ if (can_precipitate)
+ {
+ m_equations.add_equation(dof_mineral(m), &neq);
+ }
}
}
m_equations.nb_tot_variables = neq;
m_equations.nb_complementarity_variables = m_equations.nb_tot_variables - m_equations.nb_free_variables;
// Minerals => no precipitation
//
// /!\ must be set after equations ids are set
- if (not constraints.mineral_constraints.empty())
+ if (not constraints.immobile_disabled and not constraints.mineral_constraints.empty())
{
for (const MineralConstraint& v: constraints.mineral_constraints)
{
if (v.equation_type == MineralConstraintType::NoPrecipitation)
{
index_t idm = v.id_mineral;
if (idm == no_species or ideq_min(idm) == no_equation) continue;
if (not m_equations.is_box_constrained) {
m_equations.activate_box_constrained(neq);
}
m_equations.set_box_constrained(ideq_min(idm), v.param);
}
}
}
}
void AdimensionalSystem::number_eq_aqueous_component(
const AdimensionalSystemConstraints& constraints,
index_t& neq
)
{
using EqT = AqueousComponentEquationType;
// First set the charge keeper
if (constraints.charge_keeper != no_species)
{
if (constraints.charge_keeper == 0 or constraints.charge_keeper > m_data->nb_component())
{
throw std::invalid_argument("The charge keeper must be an aqueous component. Invalid argument : "
+ std::to_string(constraints.charge_keeper));
}
m_equations.type_equation(dof_component(constraints.charge_keeper)) = static_cast<index_t>(EqT::ChargeBalance);
}
// Then go over fix fugacity gas
for (const auto& it: constraints.fixed_fugacity_cs)
{
if (m_equations.type_equation(dof_component(it.id_component)) != static_cast<index_t>(EqT::NoEquation))
{
throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component)
+ "' is already constrained, a fixed fugacity condition can not be applied");
}
m_equations.type_equation(dof_component(it.id_component)) = static_cast<int>(EqT::FixedFugacity);
m_fixed_values(it.id_component) = it.log_value;
m_equations.related_species(it.id_component) = it.id_gas;
}
// Then over the fixed activity species
for (const auto& it: constraints.fixed_activity_cs)
{
if (m_equations.type_equation(dof_component(it.id_component)) != static_cast<index_t>(EqT::NoEquation))
{
throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component)
+ "' is already constrained, a fixed activity condition can not be applied.");
}
m_equations.type_equation(dof_component(it.id_component)) = static_cast<index_t>(EqT::FixedActivity);
m_fixed_values(it.id_component) = it.log_value;
}
// Then the fixed molality components
for (const auto& it: constraints.fixed_molality_cs)
{
if (m_equations.type_equation(dof_component(it.id_component)) != static_cast<index_t>(EqT::NoEquation))
{
throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component)
+ "' is already constrained, a fixed molality condition can not be applied.");
}
m_equations.type_equation(dof_component(it.id_component)) = static_cast<index_t>(EqT::FixedMolality);
m_fixed_values(it.id_component) = it.log_value;
}
// Finally number the equations
for (index_t component: m_data->range_aqueous_component())
{
// If no equation is assigned yet
if (m_equations.type_equation(dof_component(component)) == static_cast<index_t>(EqT::NoEquation))
{
// Mass is conserved for this component
//###FIXME: H[+], HO[-]
const scalar_t& total_concentration = constraints.total_concentrations(dof_component(component));
if (std::abs(total_concentration) > get_options().cutoff_total_concentration)
{
m_equations.type_equation(dof_component(component)) = static_cast<index_t>(EqT::MassConservation);
m_fixed_values(dof_component(component)) = total_concentration;
m_equations.add_equation(component, &neq);
}
else // add component to the nonactive component list
{
m_equations.add_non_active_component(component);
}
}
// If equation is already assigned
else
{
m_equations.add_equation(component, &neq);
}
}
if (stdlog::ReportLevel() >= logger::Debug and m_equations.nonactive_component.size() > 0)
{
// if in debug mode list the non active components
DEBUG << "Non active components :";
for (auto it: m_equations.nonactive_component)
{
DEBUG << " - " << it;
}
}
}
// Units
// =================
void AdimensionalSystem::set_units(const units::UnitsSet& unit_set)
{
units::UnitBaseClass::set_units(unit_set);
compute_scaling_factors();
}
void AdimensionalSystem::compute_scaling_factors()
{
// /!\ scaling factors must be set at least once per call to this functions !
// the set should be = not *= because the scaling factors are typically set twice...
// One at creation, the other before solving
// (latter is necessary because units are always reset by solver)
m_scaling_molar_volume = m_data->scaling_molar_volume(get_units());
// Unit scaling for the gaseous total concentration
// transform mol/m^3 into the right unit (pressure are always in Pa)
switch (get_units().length) {
case units::LengthUnit::decimeter:
m_scaling_gas = 1e-3;
break;
case units::LengthUnit::centimeter:
m_scaling_gas = 1e-6;
break;
default:
m_scaling_gas = 1.0;
break;
}
if (get_units().quantity == units::QuantityUnit::millimoles) {
m_scaling_gas *= 1e3;
m_scaling_molality = 1e3;
}
}
// Sums
// ==========
// ----
scalar_t AdimensionalSystem::weigthed_sum_aqueous(index_t component) const
{
scalar_t sum = 0.0;
for (index_t aqueous: m_data->range_aqueous())
{
if (not is_aqueous_active(aqueous)) continue;
sum += m_data->nu_aqueous(aqueous,component)*secondary_molality(aqueous);
}
return sum;
}
scalar_t AdimensionalSystem::diff_weigthed_sum_aqueous(index_t diff_component, index_t component) const
{
scalar_t sum = 0.0;
for (index_t aqueous: m_data->range_aqueous())
{
if (not is_aqueous_active(aqueous)) continue;
sum += m_data->nu_aqueous(aqueous,diff_component)*m_data->nu_aqueous(aqueous,component)*secondary_molality(aqueous);
}
return sum;
}
// ----
scalar_t AdimensionalSystem::weigthed_sum_sorbed(index_t component) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_sorbed(s, component)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::weigthed_sum_charge_sorbed() const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->charge_sorbed(s)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::weigthed_sum_sorbed_ssite(index_t ssite) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_ssites(s, ssite)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::diff_weigthed_sum_sorbed(index_t diff_component, index_t component) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_sorbed(s, diff_component)*m_data->nu_sorbed(s, component)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::diff_surface_weigthed_sum_sorbed(index_t diff_ssite, index_t component) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_sorbed(s, component)*m_data->nu_ssites(s, diff_ssite)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::diff_weigthed_sum_sorbed_ssite(index_t diff_component, index_t ssite) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_sorbed(s, diff_component)*m_data->nu_ssites(s, ssite)*sorbed_species_concentration(s);
}
return sum;
}
scalar_t AdimensionalSystem::diff_surface_weigthed_sum_sorbed_ssite(index_t diff_ssite, index_t ssite) const
{
scalar_t sum = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)) continue;
sum += m_data->nu_ssites(s, diff_ssite)*m_data->nu_ssites(s, ssite)*sorbed_species_concentration(s);
}
return sum;
}
// ----
scalar_t AdimensionalSystem::weigthed_sum_mineral(const Vector& x, index_t component) const
{
scalar_t sum = 0.0;
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation or m_data->nu_mineral(m, component) == 0.0) continue;
const auto concentration = volume_fraction_mineral(x, m)/molar_volume_mineral(m);
sum += m_data->nu_mineral(m, component)*concentration;
}
return sum;
}
// ----
scalar_t AdimensionalSystem::weigthed_sum_gas(index_t component) const
{
scalar_t sum = 0.0;
for (index_t k: m_data->range_gas())
{
if (not is_active_gas(k) or m_data->nu_gas(k, component) == 0.0) continue;
sum += m_data->nu_gas(k, component)*gas_concentration(k);
}
return sum;
}
scalar_t AdimensionalSystem::diff_weigthed_sum_gas(index_t diff_component, index_t component) const
{
scalar_t sum = 0.0;
for (index_t k: m_data->range_gas())
{
if (not is_active_gas(k) or m_data->nu_gas(k, component) == 0.0) continue;
sum += m_data->nu_gas(k, diff_component)*m_data->nu_gas(k, component)*gas_concentration(k);
}
return sum;
}
scalar_t AdimensionalSystem::edl_sqrt_surface_potential() const
{
scalar_t B=0;
const scalar_t RT = constants::gas_constant*temperature();
if (ionic_strength() != 0.0)
{
B = std::sqrt( 8*RT
* constants::vacuum_permittivity*constants::water_dielectric_constant_25
* ionic_strength()*density_water_SI()
);
} else {
B = std::sqrt( 8*RT
* constants::vacuum_permittivity*constants::water_dielectric_constant_25
* 0.00001*density_water_SI()
);
}
return B;
}
// ================ //
// //
// Residuals //
// //
// ================ //
scalar_t AdimensionalSystem::residual_water_conservation(const Vector& x) const
{
specmicp_assert(water_equation_type() == WaterEquationType::MassConservation);
scalar_t res = total_concentration_bc(0);
const scalar_t conc_w = density_water() * volume_fraction_water(x);
scalar_t aqconc = 1.0/m_data->molar_mass_basis(0) + weigthed_sum_aqueous(0);
if (surface_model() != SurfaceEquationType::NoEquation)
aqconc += weigthed_sum_sorbed(0);
res -= m_scaling_molality * conc_w * aqconc;
res -= weigthed_sum_mineral(x, 0);
if (m_data->nb_gas() > 0)
res -= weigthed_sum_gas(0);
if (m_equations.use_water_pressure_model and m_equations.solve_pressure_model)
{
// water pressure
const scalar_t aporosity = m_second.porosity;
scalar_t sat_w = volume_fraction_water(x) / aporosity;
if (sat_w < 1)
{
const scalar_t pressure = m_equations.water_pressure_model(sat_w);
res -= m_scaling_gas*(aporosity - volume_fraction_water(x))*(
pressure / (constants::gas_constant*temperature()));
}
}
res /= total_concentration_bc(0);
return res;
}
scalar_t AdimensionalSystem::residual_water_saturation(const Vector& x) const
{
specmicp_assert(water_equation_type() == WaterEquationType::SaturatedSystem);
scalar_t res = 1 - volume_fraction_water(x) - m_inert_volume_fraction;
for (index_t mineral: m_data->range_mineral())
{
res -= volume_fraction_mineral(x, mineral);
}
return res;
}
scalar_t AdimensionalSystem::residual_water_fixed_saturation(const Vector& x) const
{
specmicp_assert(water_equation_type() == WaterEquationType::FixedSaturation);
scalar_t porosity = 1 - m_inert_volume_fraction;
for (index_t mineral: m_data->range_mineral())
{
porosity -= volume_fraction_mineral(x, mineral);
}
scalar_t res = fixed_saturation_bc()*porosity - volume_fraction_water(x);
return res;
}
scalar_t AdimensionalSystem::residual_component(const Vector &x, index_t component) const
{
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::MassConservation);
const scalar_t conc_w = density_water()*volume_fraction_water(x);
scalar_t res = total_concentration_bc(component);
scalar_t aqconc = component_molality(x, component)
+ weigthed_sum_aqueous(component);
if (surface_model() != SurfaceEquationType::NoEquation) {
aqconc += weigthed_sum_sorbed(component);
}
res -= m_scaling_molality * conc_w * aqconc;
res -= weigthed_sum_mineral(x, component);
if (m_data->nb_gas() > 0)
res -= weigthed_sum_gas(component);
res /= total_concentration_bc(component);
return res;
}
scalar_t AdimensionalSystem::residual_fixed_activity(const Vector& x, index_t component) const
{
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedActivity);
scalar_t res = (fixed_activity_bc(component)
- log_gamma_component(component)
- log_component_molality(x, component)
);
res /= fixed_activity_bc(component);
return res;
}
scalar_t AdimensionalSystem::residual_fixed_molality(const Vector& x, index_t component) const
{
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedMolality);
scalar_t res = (fixed_molality_bc(component)
- log_component_molality(x, component)
);
res /= fixed_molality_bc(component);
return res;
}
scalar_t AdimensionalSystem::residual_fixed_fugacity(const Vector& x, index_t component) const
{
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedFugacity);
index_t id_g = m_equations.fixed_activity_species[component];
scalar_t res = fixed_fugacity_bc(component) + m_data->logk_gas(id_g);
for (index_t component: m_data->range_aqueous_component())
{
if (m_data->nu_gas(id_g, component) == 0) continue;
res -= m_data->nu_gas(id_g, component)*(
log_gamma_component(component)
+ log_component_molality(x, component)
);
}
res /= fixed_fugacity_bc(component);
return res;
}
scalar_t AdimensionalSystem::residual_mineral(const Vector& x, index_t m) const
{
specmicp_assert(ideq_min(m) != no_equation);
scalar_t res = m_data->logk_mineral(m);
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->nu_mineral(m, i) != 0)
{
const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i);
res -= m_data->nu_mineral(m, i)*log_activity_i;
}
}
if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m))
res -= m_data->nu_mineral(m, m_data->electron_index())*log_activity_electron(x);
return res;
}
scalar_t AdimensionalSystem::residual_charge_conservation(const Vector& x) const
{
scalar_t res = 0.0;
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->charge_component(i) != 0 and ideq_paq(i) != no_equation)
res += m_data->charge_component(i)*component_molality(x, i);
}
for (index_t j: m_data->range_aqueous())
{
if (m_data->charge_aqueous(j) == 0 and not is_aqueous_active(j)) continue;
res += m_data->charge_aqueous(j)*secondary_molality(j);
}
if (surface_model() != SurfaceEquationType::NoEquation)
{
res += weigthed_sum_charge_sorbed();
}
return m_scaling_molality*res;
}
scalar_t AdimensionalSystem::residual_surface(const Vector &x, index_t q) const
{
specmicp_assert(surface_model() != SurfaceEquationType::NoEquation);
const scalar_t conc_w = density_water()*volume_fraction_water(x);
const scalar_t tot_conc_surf = surface_total_concentration(q);
scalar_t res = tot_conc_surf;
res -= m_scaling_molality*conc_w*(free_sorption_site_concentration(x, q) + weigthed_sum_sorbed_ssite(q));
return res/tot_conc_surf;
}
scalar_t AdimensionalSystem::residual_surface_potential(const Vector& x) const
{
scalar_t A = density_water_SI()*volume_fraction_water(x)*weigthed_sum_charge_sorbed()
*constants::faraday_constant/sorption_surface_area();
scalar_t B = edl_sqrt_surface_potential();
const scalar_t RT = constants::gas_constant*temperature();
B *= std::sinh(surface_potential(x) * constants::faraday_constant/(2*RT));
//std::cout << (A-B) << "\n";
return (A - B)/density_water_SI();// /(density_water_SI()*volume_fraction_water(x));
}
scalar_t AdimensionalSystem::residual_electron(const Vector &x) const
{
specmicp_assert(electron_equation_type() == ElectronEquationType::Equilibrium);
const scalar_t conc_w = density_water()*volume_fraction_water(x);
scalar_t res = 0.0;
res -= conc_w * weigthed_sum_aqueous(m_data->electron_index());
if (surface_model() != SurfaceEquationType::NoEquation)
res -= conc_w * weigthed_sum_sorbed(m_data->electron_index());
res -= weigthed_sum_mineral(x, m_data->electron_index());
if (m_data->nb_gas() > 0)
res -= weigthed_sum_gas(m_data->electron_index());
return res/get_options().scaling_electron;
}
void AdimensionalSystem::get_residuals(const Vector& x, Vector& residual)
{
residual.resize(total_variables());
// initialisation of 'main' secondary variables
// They are especially nessary if the finite difference jacobian is used
// and for the linesearch
m_second.porosity = porosity(x);
set_volume_fraction_gas_phase(x);
set_secondary_concentration(x);
set_sorbed_concentrations(x);
set_pressure_fugacity(x);
//
// water
if (ideq_w() != no_equation)
{
switch (water_equation_type())
{
case WaterEquationType::MassConservation:
residual(ideq_w()) = residual_water_conservation(x);
break;
case WaterEquationType::SaturatedSystem:
residual(ideq_w()) = residual_water_saturation(x);
break;
case WaterEquationType::FixedSaturation:
residual(ideq_w()) = residual_water_fixed_saturation(x);
break;
case WaterEquationType::NoEquation:
break;
}
}
// aqueous component
for (index_t i: m_data->range_aqueous_component())
{
switch (aqueous_component_equation_type(i))
{
case AqueousComponentEquationType::NoEquation:
break;
case AqueousComponentEquationType::MassConservation:
residual(ideq_paq(i)) = residual_component(x, i);
break;
case AqueousComponentEquationType::ChargeBalance:
residual(ideq_paq(i)) = residual_charge_conservation(x);
break;
case AqueousComponentEquationType::FixedActivity:
residual(ideq_paq(i)) = residual_fixed_activity(x, i);
break;
case AqueousComponentEquationType::FixedFugacity:
residual(ideq_paq(i)) = residual_fixed_fugacity(x, i);
break;
case AqueousComponentEquationType::FixedMolality:
residual(ideq_paq(i)) = residual_fixed_molality(x, i);
break;
}
}
// surface
if (surface_model() != SurfaceEquationType::NoEquation)
{
for (index_t q: m_data->range_ssites()) {
if (ideq_surf(q) != no_equation)
residual(ideq_surf(q)) = residual_surface(x, q);
}
if (ideq_surf_pot() != no_equation) {
residual(ideq_surf_pot()) = residual_surface_potential(x);
}
}
// mineral
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m);
}
// electron
if (ideq_electron() != no_equation ) residual(ideq_electron()) = residual_electron(x);
}
// ================ //
// //
// Jacobian //
// //
// ================ //
void AdimensionalSystem::get_jacobian(Vector& x, Matrix& jacobian)
// //non-optimized Finite difference, for test only
#ifdef SPECMICP_DEBUG_EQUATION_FD_JACOBIAN
{
finite_difference_jacobian(x, jacobian);
return;
}
#else // analytical jacobian
{
analytical_jacobian(x, jacobian);
return;
}
#endif
void AdimensionalSystem::finite_difference_jacobian(Vector& x, Matrix& jacobian)
{
const int neq = total_variables();
Eigen::VectorXd res(neq);
Eigen::VectorXd perturbed_res(neq);
jacobian.setZero(neq, neq);
get_residuals(x, res);
for (int j=0; j<neq; ++j)
{
double h = 1e-8*std::abs(x(j));
if (h<1e-16) h = 1e-8;
double tmp = x(j);
x(j) += h;
h = x(j) - tmp;
get_residuals(x, perturbed_res);
for (int i=0; i<neq; ++i)
{
jacobian(i, j) = (perturbed_res(i) - res(i))/h;
}
x(j) = tmp;
}
}
void AdimensionalSystem::analytical_jacobian(Vector& x, Matrix& jacobian)
{
const int neq = total_variables();
jacobian.setZero(neq, neq);
// water
jacobian_water(x, jacobian);
// aqueous component
jacobian_aqueous_components(x, jacobian);
// surface
if (surface_model() != SurfaceEquationType::NoEquation) {
jacobian_surface(x, jacobian);
if (ideq_surf_pot() != no_equation) jacobian_surface_potential(x, jacobian);
}
// mineral equilibrium
jacobian_minerals(x, jacobian);
// electron
if (ideq_electron() != no_equation) jacobian_electron(x, jacobian);
}
void AdimensionalSystem::jacobian_water(Vector& x, Matrix& jacobian)
{
if (water_equation_type() == WaterEquationType::MassConservation)
{
const index_t idw = ideq_w();
const scalar_t rho_w = density_water();
const scalar_t factor = total_concentration_bc(0);
scalar_t tmp = - 1.0
/m_data->molar_mass_basis(0);
tmp -= weigthed_sum_aqueous(0);
tmp -= weigthed_sum_sorbed(0);
tmp *= m_scaling_molality*rho_w;
tmp -= -weigthed_sum_gas(0);
if (m_equations.use_water_pressure_model and m_equations.solve_pressure_model )
{
// The jacobian of the pressure model is
// computed using finite difference
const scalar_t sat_w = volume_fraction_water(x)/ porosity(x);
if (sat_w >= 1.0)
{
Vector residual;
get_residuals(x, residual);
// doesn't make sense to have a saturation
// bigger than one in this case
ERROR << "Saturation greater that one detected : "
<< sat_w
<< " - porosity : " << porosity(x)
<< " - tot vol frac : " << sum_volume_fraction_minerals(x)
<< ", skip pressure model \n current solution \n ---- \n"
<< x
<< "\n ---- \n"
<< "residuals \n ----- \n"
<< residual
<< "\n ------ \n";
}
else
{
// compute the perturbation
scalar_t sp = sat_w*(1.0 + eps_jacobian);
if (sp == 0.0) sp = eps_jacobian;
const scalar_t h = sp - sat_w;
// can we save one call here ? cache
const scalar_t pv_sds = m_equations.water_pressure_model(sp);
const scalar_t pv_s = m_equations.water_pressure_model(sat_w);
const scalar_t diff = (pv_sds - pv_s) / h;
// add the contribution
tmp -= m_scaling_gas/(constants::gas_constant*temperature()) *(
(1.0-sat_w)*diff - pv_s
);
}
}
jacobian(idw, idw) = tmp/factor;
// aqueous component
const scalar_t conc_w = density_water()*volume_fraction_water(x);
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation) continue;
scalar_t tmp = 0.0;
tmp -= diff_weigthed_sum_aqueous(k, 0);
tmp -= diff_weigthed_sum_sorbed(k, 0);
// fixme gas
tmp *= m_scaling_molality*conc_w;
tmp -= diff_weigthed_sum_gas(k, 0);
jacobian(idw, ideq_paq(k)) = log10*tmp/factor;
}
// electron
if (ideq_electron() != no_equation)
{
scalar_t tmp = 0.0;
tmp -= diff_weigthed_sum_aqueous(m_data->electron_index(), 0);
tmp -= diff_weigthed_sum_sorbed(m_data->electron_index(), 0);
tmp*= m_scaling_molality*conc_w;
tmp -= diff_weigthed_sum_gas(m_data->electron_index(), 0);
jacobian(idw, ideq_electron()) = log10*tmp/factor;
}
// mineral
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(idw, ideq_min(m)) = -m_data->nu_mineral(m, 0)/molar_volume_mineral(m)/factor;
}
if (surface_model() != SurfaceEquationType::NoEquation)
{
// sorption sites
for (index_t q: m_data->range_ssites())
{
if (ideq_surf(q) == no_equation) continue;
scalar_t tmp=diff_surface_weigthed_sum_sorbed(q,0);
tmp *= -m_scaling_molality*conc_w;
jacobian(idw,ideq_surf(q)) = log10*tmp/factor;
}
// surface potential
if (ideq_surf_pot() != no_equation)
{
scalar_t tmp = 0;
for (index_t p: m_data->range_sorbed()) {
if (m_data->nu_sorbed(p, 0) == 0 )continue;
tmp += m_data->nu_sorbed(p, 0) *
m_data->charge_sorbed(p) *
m_second.sorbed_concentrations(p);
}
tmp *= constants::faraday_constant/(constants::gas_constant*temperature());
tmp *= m_scaling_molality*conc_w;
jacobian(idw, ideq_surf_pot()) = tmp/factor;
}
}
}
else if (water_equation_type() == WaterEquationType::SaturatedSystem)
{
const index_t idw = ideq_w();
jacobian(idw, idw) = -1;
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(idw, ideq_min(m)) = -1;
}
}
else if (water_equation_type() == WaterEquationType::FixedSaturation)
{
const index_t idw = ideq_w();
jacobian(idw, idw) = -1;
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(idw, ideq_min(m)) = - fixed_saturation_bc();
}
}
}
void AdimensionalSystem::jacobian_aqueous_components(Vector& x, Matrix& jacobian)
{
for (index_t i: m_data->range_aqueous_component())
{
const index_t idp = ideq_paq(i);
if (idp == no_equation) continue;
switch (aqueous_component_equation_type(i))
{
case AqueousComponentEquationType::NoEquation:
continue;
// Mass balance equation
// =====================
case AqueousComponentEquationType::MassConservation:
{
const scalar_t conc_w = density_water()*volume_fraction_water(x);
const scalar_t factor = total_concentration_bc(i);
// Aqueous components
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation) continue;
scalar_t tmp_iip = 0;
if (k == i) tmp_iip -= component_molality(x, i);
tmp_iip -= diff_weigthed_sum_aqueous(k, i) + diff_weigthed_sum_sorbed(k, i);
tmp_iip *= m_scaling_molality*conc_w;
tmp_iip -= diff_weigthed_sum_gas(k, i);
jacobian(idp, ideq_paq(k)) = log10*tmp_iip/factor;
}
// Minerals
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(idp, ideq_min(m)) = - m_data->nu_mineral(m, i)/molar_volume_mineral(m)/factor;
}
// Water
if (ideq_w() != no_equation)
{
scalar_t tmp_iw = (
- component_molality(x, i)
- weigthed_sum_aqueous(i)
- weigthed_sum_sorbed(i)
);
tmp_iw *= m_scaling_molality*density_water();
jacobian(idp, ideq_w()) = tmp_iw/factor;
}
// Surface
if (surface_model() != SurfaceEquationType::NoEquation)
{
// sorption sites
for (index_t q: m_data->range_ssites())
{
if (ideq_surf(q) == no_equation) continue;
scalar_t tmp=diff_surface_weigthed_sum_sorbed(q,i);
tmp *= -m_scaling_molality*conc_w;
jacobian(idp, ideq_surf(q)) = log10*tmp/factor;
}
// surface potential
if (ideq_surf_pot() != no_equation)
{
scalar_t tmp = 0.0;
for (index_t p: m_data->range_sorbed()) {
if (not is_active_sorbed(p) or m_data->nu_sorbed(p, i) == 0) continue;
tmp += -m_data->nu_sorbed(p, i) *
m_data->charge_sorbed(p) *
m_second.sorbed_concentrations(p);
}
tmp *= constants::faraday_constant/(
constants::gas_constant*temperature());
tmp *= -m_scaling_molality*conc_w;
jacobian(idp, ideq_surf_pot()) = tmp/factor;
}
}
if (ideq_electron() != no_equation)
{
scalar_t tmp = 0.0;
tmp -= diff_weigthed_sum_aqueous(m_data->electron_index(), i);
tmp -= diff_weigthed_sum_sorbed(m_data->electron_index(), i);
tmp*= m_scaling_molality*conc_w;
tmp -= diff_weigthed_sum_gas(m_data->electron_index(), i);
jacobian(idp, ideq_electron()) = log10*tmp/factor;
}
break;
}
// Charge balance equation
// =======================
case AqueousComponentEquationType::ChargeBalance:
{
// Aqueous components
for (index_t k: m_data->range_aqueous_component())
{
const index_t idc = ideq_paq(k);
if (idc == no_equation) continue;
scalar_t tmp_drdb = 0.0;
if (m_data->charge_component(k) != 0.0)
tmp_drdb = m_data->charge_component(k)*component_molality(x,k);
// Secondary species
for (index_t j: m_data->range_aqueous())
{
if ( not is_aqueous_active(j)
or m_data->nu_aqueous(j, k) == 0.0
or m_data->charge_aqueous(j) == 0.0
)
continue;
tmp_drdb += m_data->nu_aqueous(j, k)
* m_data->charge_aqueous(j)
* secondary_molality(j);
}
// sorbed species
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)
or m_data->nu_sorbed(s, k) == 0.0
or m_data->charge_sorbed(s) == 0.0
) continue;
tmp_drdb += m_data->nu_sorbed(s, k)
* m_data->charge_sorbed(s)
* sorbed_species_concentration(s);
}
jacobian(idp, idc) = m_scaling_molality*log10*tmp_drdb;
}
if (surface_model() != SurfaceEquationType::NoEquation)
{
// Sorption sites
for (index_t q: m_data->range_ssites())
{
const index_t idq = ideq_surf(q);
if (idq == no_equation) continue;
scalar_t tmp_drdlsq = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)
or m_data->nu_ssites(s, q) == 0.0
or m_data->charge_sorbed(s) == 0.0
) continue;
tmp_drdlsq += m_data->nu_ssites(s, q)
* m_data->charge_sorbed(s)
* sorbed_species_concentration(s);
}
jacobian(idp, idq) = m_scaling_molality*log10*tmp_drdlsq;
}
// Surface potential
if (ideq_surf_pot() != no_equation) {
scalar_t tmp_drdpsi = 0.0;
for (index_t s: m_data->range_sorbed())
{
if (not is_active_sorbed(s)
or m_data->charge_sorbed(s) == 0.0
) continue;
tmp_drdpsi += - std::pow(m_data->charge_sorbed(s), 2)
* sorbed_species_concentration(s);
}
jacobian(idp, ideq_surf_pot()) = m_scaling_molality*tmp_drdpsi*
constants::faraday_constant /
(constants::gas_constant*temperature());
}
}
break;
}
// Fixed activity equation
// =======================
case AqueousComponentEquationType::FixedActivity:
{
jacobian(idp, idp) = -1.0/fixed_activity_bc(i);
break;
}
// Fixed fugacity equation
// =======================
case AqueousComponentEquationType::FixedFugacity:
{
index_t id_g = m_equations.fixed_activity_species[i];
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation or m_data->nu_gas(id_g, k) == 0.0) continue;
jacobian(idp, ideq_paq(k)) = -m_data->nu_gas(id_g, k)/fixed_fugacity_bc(i);
}
break;
} // end case
// Fixed molality component
// ========================
case AqueousComponentEquationType::FixedMolality:
{
jacobian(idp, idp) = -1.0/fixed_molality_bc(i);
break;
}
} // end switch
}
}
void AdimensionalSystem::jacobian_minerals(Vector& x, Matrix& jacobian)
{
for (index_t m: m_data->range_mineral())
{
const index_t idm = ideq_min(m);
if (idm == no_equation) continue;
for (index_t i: m_data->range_aqueous_component())
{
if (ideq_paq(i) == no_equation) continue;
jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i);
}
if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m))
jacobian(idm, ideq_electron()) = -m_data->nu_mineral(m, m_data->electron_index());
}
}
void AdimensionalSystem::jacobian_surface(Vector& x, Matrix& jacobian)
{
for (index_t q: m_data->range_ssites())
{
const index_t idq = ideq_surf(q);
if (idq == no_equation) continue;
const scalar_t factor = m_scaling_molality*density_water()/surface_total_concentration(q);
// Water
if (ideq_w() != no_equation) {
scalar_t tmp = free_sorption_site_concentration(x, q)+weigthed_sum_sorbed_ssite(q);
jacobian(idq, ideq_w()) = -factor*tmp;
}
// Aqueous component
for (index_t i: m_data->range_aqueous_component()) {
const index_t idc = ideq_paq(i);
if (idc == no_equation) continue;
scalar_t tmp = log10*volume_fraction_water(x)*diff_weigthed_sum_sorbed_ssite(i, q);
jacobian(idq, idc) = -factor*tmp;
}
// Surface site
for (index_t qp: m_data->range_ssites()) {
const index_t idqp = ideq_surf(qp);
if (idqp == no_equation) continue;
scalar_t tmp = diff_surface_weigthed_sum_sorbed_ssite(qp, q);
if (qp == q) {
tmp += free_sorption_site_concentration(x, q);
}
jacobian(idq, idqp) = -factor*volume_fraction_water(x)*log10*tmp;
}
// Surface potential
if (ideq_surf_pot() != no_equation) {
scalar_t tmp = 0;
for (index_t p: m_data->range_sorbed()) {
if (not is_active_sorbed(p) or m_data->nu_ssites(p, q) == 0) continue;
tmp += m_data->nu_ssites(p, q) *
m_data->charge_sorbed(p) *
m_second.sorbed_concentrations(p);
}
tmp *= - constants::faraday_constant/(constants::gas_constant*temperature());
jacobian(idq, ideq_surf_pot()) = -factor*volume_fraction_water(x)*tmp;
}
// Mineral
// Electron
}
}
void AdimensionalSystem::jacobian_surface_potential(Vector& x, Matrix& jacobian)
{
const index_t idqs = ideq_surf_pot();
const scalar_t factor_A = constants::faraday_constant/sorption_surface_area();
specmicp_assert(idqs != no_equation);
// water
if (ideq_w() != no_equation) {
jacobian(idqs, ideq_w()) = factor_A*weigthed_sum_charge_sorbed();
}
// aq. compoments
for (index_t i: m_data->range_aqueous_component()) {
const index_t idc = ideq_paq(i);
if (idc == no_equation) continue;
scalar_t tmp = 0;
for (index_t p: m_data->range_sorbed()) {
if (m_data->nu_sorbed(p, i) == 0 ) continue;
tmp += m_data->nu_sorbed(p, i) *
m_data->charge_sorbed(p) *
m_second.sorbed_concentrations(p);
}
jacobian(idqs, idc) = factor_A*volume_fraction_water(x)*log10*tmp;
}
// surface site
{
for (index_t q: m_data->range_ssites()) {
if (ideq_surf(q) == no_equation) continue;
scalar_t tmp = 0.0;
for (index_t p: m_data->range_sorbed()) {
if (not is_active_sorbed(p) or m_data->nu_ssites(p, q)==0.0) continue;
tmp += m_data->charge_sorbed(p) * m_data->nu_ssites(p, q)*m_second.sorbed_concentrations(p);
}
jacobian(idqs, ideq_surf(q)) = factor_A*volume_fraction_water(x)*log10*tmp;
}
}
// surface potential
{
const scalar_t RT = constants::gas_constant*temperature();
scalar_t dB = constants::faraday_constant/(2*RT)*edl_sqrt_surface_potential();
dB *= std::cosh(surface_potential(x)*constants::faraday_constant/(2*RT) );
//
scalar_t dA = 0.0;
for (index_t p: m_data->range_sorbed()) {
if (not is_active_sorbed(p)) continue;
dA += -std::pow(m_data->charge_sorbed(p),2) *
m_second.sorbed_concentrations(p);
}
dA *= factor_A*density_water_SI()*volume_fraction_water(x)*constants::faraday_constant/(RT);
jacobian(idqs, idqs) = (dA-dB)/(density_water_SI());
}
}
void AdimensionalSystem::jacobian_electron(Vector& x, Matrix& jacobian)
{
const auto ide = ideq_electron();
const auto dofe = m_data->electron_index();
const scalar_t conc_w = density_water()*volume_fraction_water(x);
const scalar_t factor = get_options().scaling_electron;
// Aqueous components
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation) continue;
scalar_t tmp_eip = 0;
tmp_eip -= diff_weigthed_sum_aqueous(k, dofe);
tmp_eip -= diff_weigthed_sum_sorbed(k, dofe);
tmp_eip *= conc_w;
tmp_eip -= diff_weigthed_sum_gas(k, dofe);
jacobian(ide, ideq_paq(k)) = log10*tmp_eip/factor;
}
// Minerals
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(ide, ideq_min(m)) = - m_data->nu_mineral(m, dofe)/molar_volume_mineral(m)/factor;
}
// Water
if (ideq_w() != no_equation)
{
scalar_t tmp_iw = 0;
tmp_iw -= weigthed_sum_aqueous(dofe);
tmp_iw -= weigthed_sum_sorbed(dofe);
tmp_iw *= density_water();
jacobian(ide, ideq_w()) = tmp_iw/factor;
}
// Surface
// if (ideq_surf() != no_equation)
// {
// scalar_t tmp_s = -conc_w*diff_surface_weigthed_sum_sorbed(dofe);
// jacobian(ide, ideq_surf()) = tmp_s/factor;
// }
// Electron
if (ideq_electron() != no_equation)
{
scalar_t tmp = 0.0;
tmp -= diff_weigthed_sum_aqueous(dofe, dofe);
tmp -= diff_weigthed_sum_sorbed(dofe, dofe);
tmp*= conc_w;
tmp -= diff_weigthed_sum_gas(dofe, dofe);
jacobian(ide, ide) = log10*tmp/factor;
}
}
// ========================== //
// //
// Secondary variables //
// //
// ========================== //
void AdimensionalSystem::set_secondary_variables(const Vector& x)
{
m_second.porosity = 1 - sum_volume_fraction_minerals(x) - m_inert_volume_fraction;
if (surface_model() != SurfaceEquationType::NoEquation) {
set_sorbed_concentrations(x);
}
set_secondary_concentration(x);
compute_log_gamma(x);
set_volume_fraction_gas_phase(x);
set_pressure_fugacity(x);
}
void AdimensionalSystem::set_volume_fraction_gas_phase(const Vector& x)
{
m_second.volume_fraction_gas = m_second.porosity - volume_fraction_water(x);
}
void AdimensionalSystem::set_pressure_fugacity(const Vector& x)
{
const auto rt = constants::gas_constant*temperature();
for (index_t k: m_data->range_gas())
{
if (not is_active_gas(k)) continue;
scalar_t logp = -m_data->logk_gas(k);
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->nu_gas(k, i) == 0.0) continue;
const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i);
logp += m_data->nu_gas(k, i) * log_activity_i;
}
if (ideq_electron() != no_equation and m_data->is_gas_half_cell_reaction(k))
logp += m_data->nu_gas(k, m_data->electron_index())*log_activity_electron(x);
m_second.gas_fugacity(k) = pow10(logp);
const scalar_t pressure = gas_fugacity(k)*gas_total_pressure();
const scalar_t concentration = m_scaling_gas*volume_fraction_gas_phase()*pressure/rt;
m_second.gas_concentration(k) = concentration;
}
}
void AdimensionalSystem::set_secondary_concentration(const Vector& x)
{
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j))
{
m_second.secondary_molalities(j) = 0.0;
continue;
}
scalar_t logconc = -m_data->logk_aqueous(j) - log_gamma_secondary(j);
for (index_t k: m_data->range_aqueous_component())
{
if (m_data->nu_aqueous(j, k) == 0) continue;
const auto log_activity_k = log_component_molality(x, k) + log_gamma_component(k);
logconc += m_data->nu_aqueous(j, k)*log_activity_k;
}
if (ideq_electron() != no_equation and m_data->is_half_cell_reaction(j))
logconc += m_data->nu_aqueous(j, m_data->electron_index())*log_activity_electron(x);
m_second.secondary_molalities(j) = pow10(logconc);
}
}
void AdimensionalSystem::set_sorbed_concentrations(const Vector& x)
{
for (index_t s: m_data->range_sorbed())
{
index_t q = m_data->sorbed.surface_site(s);
if (not is_active_sorbed(s))
{
m_second.sorbed_concentrations(s) = 0.0;
continue;
}
scalar_t logconc = -m_data->logk_sorbed(s) +
m_data->sorbed.nu_ssite(s)*log_free_sorption_site_concentration(x, q);
for (index_t k: m_data->range_aqueous_component())
{
if (m_data->nu_sorbed(s, k) != 0.0) {
const auto log_activity_k = log_component_molality(x, k) + log_gamma_component(k);
logconc += m_data->nu_sorbed(s, k)*log_activity_k;
}
}
if (ideq_electron() != no_equation and m_data->is_sorbed_half_cell_reaction(s))
logconc += m_data->nu_sorbed(s, m_data->electron_index())*log_activity_electron(x);
m_second.sorbed_concentrations(s) = pow10(logconc);
if (surface_model() == SurfaceEquationType::EDL) {
m_second.sorbed_concentrations(s) *= std::exp(
-m_data->charge_sorbed(s)*surface_potential(x)*constants::faraday_constant /
(constants::gas_constant*temperature()));
}
}
}
void AdimensionalSystem::set_ionic_strength(const Vector& x)
{
scalar_t ionic = 0;
for (index_t i: m_data->range_aqueous_component())
{
if (ideq_paq(i) == no_equation or m_data->charge_component(i) == 0) continue;
ionic += component_molality(x, i)*std::pow(m_data->charge_component(i),2);
}
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j) or m_data->charge_aqueous(j) == 0) continue;
ionic += secondary_molality(j)*std::pow(m_data->charge_aqueous(j),2);
}
ionic_strength() = ionic/2;
}
void AdimensionalSystem::compute_log_gamma(const Vector& x)
{
set_ionic_strength(x);
const scalar_t ion_str = ionic_strength();
const scalar_t sqrti = std::sqrt(ionic_strength());
vector_debye_huckel_component(ion_str, sqrti);
vector_debye_huckel_aqueous(ion_str, sqrti);
}
void AdimensionalSystem::vector_debye_huckel_component(
scalar_t ionic_strength,
scalar_t sqrt_ionic
)
{
Vector& log_gamma = m_second.loggamma;
Matrix& ionic_param = m_data->components.m_ionic_param.m_matrix;
using IoParam = database::IonicModelParameters;
for (auto i: m_data->range_aqueous_component())
{
const scalar_t zisquare = std::pow(ionic_param(i, IoParam::charge_ind), 2);
log_gamma(i) = - ( constants::Adebye * zisquare * sqrt_ionic);
log_gamma(i) /= ( 1 + ionic_param(i, IoParam::adebye_ind)*constants::Bdebye*sqrt_ionic);
log_gamma(i) += ionic_param(i, IoParam::bdebye_ind)*ionic_strength;
}
}
void AdimensionalSystem::vector_debye_huckel_aqueous(
scalar_t ionic_strength,
scalar_t sqrt_ionic
)
{
Vector& log_gamma = m_second.loggamma;
Matrix& ionic_param = m_data->aqueous.m_ionic_param.m_matrix;
using IoParam = database::IonicModelParameters;
const index_t offset = m_data->nb_component();
for (auto j: m_data->range_aqueous())
{
const scalar_t zisquare = std::pow(ionic_param(j, IoParam::charge_ind), 2);
log_gamma(offset+j) = - (constants::Adebye * zisquare * sqrt_ionic);
log_gamma(offset+j) /= ( 1 + ionic_param(j, IoParam::adebye_ind)*constants::Bdebye*sqrt_ionic);
log_gamma(offset+j) += ionic_param(j, IoParam::bdebye_ind)*ionic_strength;
}
}
bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual)
{
if (not get_options().non_ideality)
{
set_secondary_variables(x);
// we still need to compute secondary species !
// if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation)
// {
// set_secondary_variables(x);
// }
return true;
}
not_in_linesearch = true;
scalar_t previous_norm = m_second.loggamma.norm();
if (previous_norm == 0) previous_norm = 1.0;
bool may_have_converged = false;
if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation)
{
m_equations.solve_pressure_model = true;
// Use fixed point iterations for non-ideality
for (int i=0; i<get_options().non_ideality_max_iter; ++i)
{
set_secondary_variables(x);
compute_log_gamma(x);
// convergence check
if ( (std::abs(previous_norm - m_second.loggamma.norm())/previous_norm <
get_options().non_ideality_tolerance)
) {
may_have_converged = true;
break;
}
previous_norm = m_second.loggamma.norm();
}
}
else {
m_equations.solve_pressure_model = false;
}
return may_have_converged;
}
double AdimensionalSystem::max_lambda(const Vector& x, const Vector& update)
{
if (ideq_w() != no_equation)
{
return 1.0/std::max(
1.0, -update(0)/(get_options().under_relaxation_factor*x(0))
);
}
else
{
return 1.0;
}
}
AdimensionalSystemSolution AdimensionalSystem::unsafe_get_solution(
Vector& xtot,
const Vector& x
)
{
return AdimensionalSystemSolution(xtot,
m_second.secondary_molalities,
m_second.loggamma,
m_second.ionic_strength,
m_second.gas_fugacity,
m_second.sorbed_concentrations,
m_inert_volume_fraction);
}
AdimensionalSystemSolution AdimensionalSystem::get_solution(
Vector& xtot,
const Vector& x
)
{
// make sure secondary variables are ok
double previous_norm = m_second.loggamma.norm();
set_volume_fraction_gas_phase(x);
set_pressure_fugacity(x);
set_secondary_concentration(x);
if (surface_model() != SurfaceEquationType::NoEquation) {
set_sorbed_concentrations(x);
}
if (get_options().non_ideality)
{
compute_log_gamma(x);
double error = std::abs(previous_norm - m_second.loggamma.norm());
if (error > 1e-6)
{
WARNING << "Activity coefficient have not converged !" << std::endl
<< "output can not be trusted\n Difference : "
+std::to_string(error);
}
}
// Set the correct value for the water volume fraction
if (ideq_w() == no_equation)
{
xtot(dof_water()) = volume_fraction_water(x);
}
return unsafe_get_solution(xtot, x);
}
// Water, saturation and density
// ==============================
scalar_t AdimensionalSystem::density_water() const {
return laws::density_water(get_units());
}
scalar_t AdimensionalSystem::density_water_SI() const {
return laws::density_water(units::SI_units);
}
scalar_t AdimensionalSystem::volume_fraction_water(const Vector& x) const
{
if (ideq_w() != no_equation)
return x(ideq_w());
else
return porosity(x);
}
scalar_t AdimensionalSystem::volume_fraction_mineral(const Vector& x, index_t mineral) const
{
specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral());
if (ideq_min(mineral) == no_equation) return 0.0;
else return x(ideq_min(mineral));
}
scalar_t AdimensionalSystem::sum_volume_fraction_minerals(const Vector& x) const
{
scalar_t sum_saturations = 0.0;
for (index_t mineral: m_data->range_mineral())
{
sum_saturations += volume_fraction_mineral(x, mineral);
}
return sum_saturations;
}
// Starting guess
// ==============
void AdimensionalSystem::reasonable_starting_guess(Vector &xtot)
{
xtot.resize(total_dofs());
xtot(dof_water()) = get_options().restart_water_volume_fraction;
for (index_t i: m_data->range_aqueous_component())
{
xtot(dof_component(i)) = -6.0;
}
if (surface_model() != SurfaceEquationType::NoEquation) {
for (index_t q: m_data->range_ssites()) {
if (ideq_surf(q) != no_equation)
xtot(dof_surface(q)) = std::log10(0.1*surface_total_concentration(q));
else
xtot(dof_surface(q)) = -HUGE_VAL;
}
if (surface_model() == SurfaceEquationType::EDL) {
xtot(dof_surface_potential()) = 1e-3;
} else {
xtot(dof_surface_potential()) = -HUGE_VAL;
}
}
if (ideq_electron() != no_equation)
xtot(dof_electron()) = -4;
else
xtot(dof_electron()) = -HUGE_VAL;
xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero();
m_second = SecondaryVariables(m_data);
}
void AdimensionalSystem::reasonable_restarting_guess(Vector& xtot)
{
static std::mt19937 gen(std::random_device{}());
std::uniform_real_distribution<> dis(-2, 2);
xtot(dof_water()) = get_options().restart_water_volume_fraction;
for (index_t i: m_data->range_aqueous_component())
{
//if (xtot(dof_component(i)) > 0 or xtot(dof_component(i)) < -9)
xtot(i) = get_options().restart_concentration + dis(gen);
}
if (surface_model() != SurfaceEquationType::NoEquation) {
for (index_t q: m_data->range_ssites()) {
if (ideq_surf(q) != no_equation)
xtot(dof_surface(q)) = std::log10(0.1*surface_total_concentration(q));
else
xtot(dof_surface(q)) = -HUGE_VAL;
}
if (surface_model() == SurfaceEquationType::EDL) {
xtot(dof_surface_potential()) = 1e-3;
} else {
xtot(dof_surface_potential()) = -HUGE_VAL;
}
}
if (ideq_electron() != no_equation)
xtot(dof_electron()) = -4;
else
xtot(dof_electron()) = -HUGE_VAL;
xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero();
m_second = SecondaryVariables(m_data);
}
} // end namespace specmicp
diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp
index e83c50c..8ada4a7 100644
--- a/src/specmicp/adimensional/adimensional_system.hpp
+++ b/src/specmicp/adimensional/adimensional_system.hpp
@@ -1,883 +1,881 @@
/* =============================================================================
-
- Copyright (c)
- 2014 - 2017 F. Georget <fabieng@princeton.edu> Princeton University
- 2017 - 2018 F. Georget <fabien.georget@epfl.ch> EPFL
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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. *
============================================================================= */
#ifndef SPECMICP_ADIMENSIONALSYSTEM_HPP
#define SPECMICP_ADIMENSIONALSYSTEM_HPP
//! \file adimensional_system.hpp The MiCP program to solve speciation
#include "adimensional_system_numbering.hpp"
#include "adimensional_system_structs.hpp"
#include "specmicp_common/types.hpp"
#ifndef SPECMICP_DATABASE_HPP
#include "specmicp_database/database.hpp"
#endif
#include "specmicp_common/micpsolver/micpprog.hpp"
#include "specmicp_common/physics/units.hpp"
#include "specmicp_common/options_handler.hpp"
//! \namespace specmicp main namespace for the SpecMiCP solver
namespace specmicp {
struct AdimensionalSystemSolution;
//! \brief The equilibrium speciation problem
//!
//! Represent the equilibrium speciation problem as a Mixed Complementarity program
//!
//! The main variables are :
//! - \f$phi_w\f$ the volume fraction of water (total saturation)
//! - \f$b_i\f$ the molality of aqueous component
//! - \f$phi_m\f$ the volume fraction of mineral
//! - \f$s_f\f$ the concentration (molality) of free sorption sites
//!
//! The solid phase equilibrium is solved by using a complementarity formulation
//! \f[ S^t_m \geq 0, \; - \mathrm{SI}_m \geq 0 , \; \mathrm{and} \; - S^t_m \mathrm{SI}_m = 0 \f]
//!
//! The secondary variables (molalities of secondary aqueous species, gas fugacities, ionic strength, ...)
//! are solved using a fixed point iteration at the beginning of each iteration.
//!
//! This class should only be used through the AdimensionalSystemSolver.
class SPECMICP_DLL_LOCAL AdimensionalSystem :
public micpsolver::MiCPProg<AdimensionalSystem>,
public AdimemsionalSystemNumbering,
public OptionsHandler<AdimensionalSystemOptions>,
private units::UnitBaseClass
{
public:
//! \brief Create a reduced system
//!
//! \param ptrdata Shared pointer to the thermodynamic database
//! \param constraints Constraints to apply to the system
//! \param options Numerical options, mainly related to the secondary variables
//! \param units_set The set of units
AdimensionalSystem(RawDatabasePtr& ptrdata,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemOptions& options=AdimensionalSystemOptions(),
const units::UnitsSet& units_set=units::UnitsSet()
);
//! \brief Create a reduced system, initialize the system using a previous solution
//!
//! \param ptrdata Shared pointer to the thermodynamic database
//! \param constraints Constraints to apply to the system
//! \param previous_solution The previous solution to use for initialization
//! \param options Numerical options, mainly related to the secondary variables
//! \param units_set The set of units
AdimensionalSystem(
RawDatabasePtr& ptrdata,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolution& previous_solution,
const AdimensionalSystemOptions& options=AdimensionalSystemOptions(),
const units::UnitsSet& units_set=units::UnitsSet()
);
// Variables
// ==========
//! \name Variables
//! \brief How many ?
//!
//! @{
//! \brief Return the total number of variables
index_t total_variables() const {return m_equations.nb_tot_variables;}
//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
index_t nb_free_variables() const {return m_equations.nb_free_variables;}
//! \brief Return the number of variables subjected to the complementarity conditions
index_t nb_complementarity_variables() const {return m_equations.nb_complementarity_variables;}
//! @}
//! \brief Return a pointer to the database
RawDatabasePtr get_database() {return m_data;}
// The linear system
// ==================
//! \name Residuals and Jacobian
//! \brief Compute the residuals and the jacobian
//!
//! @{
//! \brief Return the residuals
//!
//! \param x Vector of the main variables
//! \param[out] residual Vector containing the residuals
void get_residuals(const Vector& x, Vector& residual);
//! \brief Return the jacobian
//!
//! \param x Vector of the main variables
//! \param[out] jacobian Dense matrix representing the jacobian
void get_jacobian(Vector& x,
Matrix& jacobian);
//! \brief Return the residual for the water conservation equation
//!
//! \param x Vector of the main variables
scalar_t residual_water_conservation(const Vector& x) const;
//! \brief Return the residual for the water saturation equation
//!
//! \param x Vector of the main variables
scalar_t residual_water_saturation(const Vector& x) const;
//! \brief Return the residual for the water fixed saturation equation
//!
//! \param x Vector of the main variables
scalar_t residual_water_fixed_saturation(const Vector& x) const;
//! \brief Return the residual for the conservation equation of component (i)
//!
//! For water (i==0) use the method : 'residual_water'
//! \param x Vector of the main variables
//! \param i Index of the component (in the database)
scalar_t residual_component(const Vector& x, index_t i) const;
//! \brief Compute the residual for the surface sorption
//! \param x Vector of the main variables
scalar_t residual_surface(const Vector& x, index_t q) const;
//! \brief Compute the residual for the surface sorption potential
//! \param x Vector of the main variables
scalar_t residual_surface_potential(const Vector& x) const;
//! \brief Equilibrium condition for the minerals
//! \param x Vector of the main variables
//! \param m Index of the mineral (in the database)
scalar_t residual_mineral(const Vector& x, index_t m) const;
//! \brief Return the residual of the charge conservation equation
//!
//! This equation may be used instead of a mass conservation equation to
//! explicitely enforce the electroneutrality.
//! The component corresponding to this equation is called the charge keeper
//! \param x Vector of the main variables
scalar_t residual_charge_conservation(const Vector& x) const;
//! \brief Return the residual corresponding to a fixed fugacity equation
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t residual_fixed_fugacity(const Vector& x, index_t component) const;
//! \brief Return the residual corresponding to a fixed activity equation
//!
//! If 'component ' is charged, then a charge keeper should be set.
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t residual_fixed_activity(const Vector& x, index_t component) const;
//! \brief Return the residual corresponding to a fixed molality equation
//!
//! If 'component ' is charged, then a charge keeper should be set.
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t residual_fixed_molality(const Vector& x, index_t component) const;
//! \brief Return the residual corresponding to the electron equation
scalar_t residual_electron(const Vector &x) const;
//! \brief Compute the jacobian analytically
void analytical_jacobian(Vector& x, Matrix& jacobian);
//! \brief Compute the jacobian using finite difference
void finite_difference_jacobian(Vector& x, Matrix& jacobian);
//! \brief Return a scale factor to avoid negative mass during Newton's iteration
//!
//! \param x Vector of the main variables
//! \param update Update of the current iteration
double max_lambda(const Vector &x, const Vector &update);
//! @}
// Getters
// #######
// Equation id access
// ==================
//! \name Equation Id
//! \brief The equation id is the row of the equation in the residuals/jacobian
//!
//! The degree of freedom getter function are inherited from specmicp::AdimensionalSystemNumbering
//! @{
//! \brief Return the equation id
//! \param i Index of the main variable (in the set of the main variables)
index_t ideq(index_t i) const {return m_equations.ideq[i];}
//! \brief Return the equation number for conservation of water
index_t ideq_w() const {return m_equations.ideq[dof_water()];}
//! \brief Return the equation number of component 'i'
//! \param i Index of the component (in the database)
index_t ideq_paq(index_t i) const {
specmicp_assert(i < m_data->nb_component());
return m_equations.ideq[dof_component(i)];
}
//! \brief Return the equation number of the electron equation
index_t ideq_electron() const {
return m_equations.ideq[dof_electron()];
}
//! \brief Return the equation number of the free surface concentration
index_t ideq_surf(index_t q) const {
return m_equations.ideq[dof_surface(q)];
}
index_t ideq_surf_pot() const {
return m_equations.ideq[dof_surface_potential()];
}
//! \brief Return the equation number of mineral 'm'
//! \param m Index of the mineral (in the database)
index_t ideq_min(index_t m) const {
specmicp_assert(m < m_data->nb_mineral());
return m_equations.ideq[dof_mineral(m)];
}
//! @}
//! \name Equation type
//! \brief Which equation is actually solved
//!
//! @{
//! \brief Return the type of equation used for the solvent
//!
//! Can be no equation, mass conservation, saturation of the system, ...
//!
//! \sa #WaterEquationType, #aqueous_component_equation_type, #AqueousComponentEquationType
WaterEquationType water_equation_type() const {
return static_cast<WaterEquationType>(m_equations.type_equation(dof_water()));
}
//! \brief Return the type of equation for aqueous species 'component'
//!
//! For water used the method #water_equation_type
//! \param component Index of the aqueous component (in the database)
//!
//! \sa #AqueousComponentEquationType, water_equation_type, #WaterEquationType
AqueousComponentEquationType aqueous_component_equation_type(index_t component) const {
specmicp_assert(component > 0 and component < m_data->nb_component());
return static_cast<AqueousComponentEquationType>(m_equations.type_equation(dof_component(component)));
}
//! \brief Return the type of equation for the electron
ElectronEquationType electron_equation_type() const {
return static_cast<ElectronEquationType>(m_equations.type_equation(dof_electron()));
}
//! \brief Return true if an equation exist for this component
//!
//! \param component Index of the component (in the database)
bool is_active_component(index_t component) const {
specmicp_assert(component < m_data->nb_component() and component > no_species);
return m_equations.id_equation(dof_component(component)) != no_equation;
}
//! \brief Return the surface model
SurfaceEquationType surface_model() const {
return static_cast<SurfaceEquationType>(m_equations.type_equation(dof_surface_potential()));
}
//! \brief Return true if 'component' is the charge keeper
//!
//! The charge keeper is the component added or removed to satisfy the charge balance
//! \param component Index of the aqueous component (in the database)
bool is_charge_keeper(index_t component)
{
return (aqueous_component_equation_type(component) == AqueousComponentEquationType::ChargeBalance);
}
//! \brief Return the component that does not exist in the solution (inactive dof)
const std::vector<index_t>& get_non_active_component() {return m_equations.nonactive_component;}
//! @}
// Variables
// =========
//! \name2 Main variables
//! \brief Access to the main variables
//!
//! @{
//! \brief Return the volume fraction of water
//!
//! \param x Vector of the main variables
scalar_t volume_fraction_water(const Vector& x) const;
//! \brief Return the density of water in system units
scalar_t density_water() const;
//! \brief Return the density of water in SI units
scalar_t density_water_SI() const;
//! \brief Return the volume fraction of a mineral
//!
//! \param x Vector of the main variables
//! \param mineral Index of the mineral (in the database)
scalar_t volume_fraction_mineral(const Vector& x, index_t mineral) const;
//! \brief Return the volume of minerals
//!
//! \param mineral Index of the mineral (in the database)
scalar_t molar_volume_mineral(index_t mineral) const {
return m_scaling_molar_volume*m_data->unsafe_molar_volume_mineral(mineral);
}
//! \brief Return the sum of saturation of the minerals
//!
//! This corresponds to the volume fraction of the solid phases (minus the inert phase)
//! \param x Vector of the main variables
scalar_t sum_volume_fraction_minerals(const Vector& x) const;
//! \brief Return true if the mineral is active
bool is_active_mineral(index_t mineral) const {
- return (ideq_min(mineral) != no_species);
+ return (ideq_min(mineral) != no_equation);
}
//! \brief Return the porosity
//!
//! Computed 'on-the-fly'
//!
//! \param x Vector of the main variables
scalar_t porosity(const Vector& x) const {
return 1-sum_volume_fraction_minerals(x)-m_inert_volume_fraction;
}
//! \brief Return the log_10 of the component molality
//!
//! \param x Vector of the main variables
//! \param component Index of the aqueous component (in the database)
scalar_t log_component_molality(const Vector& x, index_t component) const {
specmicp_assert(ideq_paq(component) != no_equation and component < m_data->nb_component());
return x(ideq_paq(component));
}
//! \brief Return the molality of 'component'
//!
//! \param x Vector of the main variables
//! \param component Index of the aqueous component (in the database)
scalar_t component_molality(const Vector& x, index_t component) const {
return pow10(log_component_molality(x, component));
}
//! \brief Return the concentration of free sorption site
scalar_t free_sorption_site_concentration(const Vector& x, index_t q) const {
return pow10(log_free_sorption_site_concentration(x, q));
}
//! \brief Return the log_10 of the free sorption site concentration
scalar_t log_free_sorption_site_concentration(const Vector& x, index_t q) const {
specmicp_assert(ideq_surf(q) != no_equation);
return x(ideq_surf(q));
}
//! \brief Return the surface potential (in V)
scalar_t surface_potential(const Vector& x) const {
specmicp_assert(ideq_surf_pot() != no_equation);
return x(ideq_surf_pot());
}
//! \brief log activity of the electron
scalar_t log_activity_electron(const Vector& x) const {
specmicp_assert(ideq_electron() != no_equation);
return x(ideq_electron());
}
//! \brief return the activity of the electro
scalar_t activity_electron(const Vector& x) const {
return pow10(log_activity_electron(x));
}
//! @}
//! \name Secondary variables
//! \brief Access to the secondary variables
//!
//! @{
//! \brief Return the concentration of secondary species 'aqueous'
//!
//! \param aqueous Index of the secondary aqueous species (in the database)
scalar_t secondary_molality(index_t aqueous) const {
if (not is_aqueous_active(aqueous)) return 0.0;
return m_second.secondary_molalities(dof_aqueous(aqueous));
}
//! \brief Return true if 'aqueous' is active
//!
//! i.e. Return true if all its component are present in the system
//! \param aqueous Index of the secondary aqueous species (in the database)
bool is_aqueous_active(index_t aqueous) const {
return m_equations.active_aqueous[dof_aqueous(aqueous)];
}
//! \brief Return log_10(γ_i) where i is a component
//!
//! γ is the activity coefficient
//! \param component Index of the aqueous component (in the database)
scalar_t log_gamma_component(index_t component) const {
return m_second.loggamma(dof_component_gamma(component));
}
//! \brief Return log_10(γ_j) where j is a secondary aqueous species
//!
//! γ is the activity coefficient
//! \param aqueous Index of the secondary aqueous species (in the database)
scalar_t log_gamma_secondary(index_t aqueous) const {
return m_second.loggamma(dof_aqueous_gamma(aqueous));
}
//! \brief Return the ionic strength
scalar_t ionic_strength() const noexcept {
return m_second.ionic_strength;
}
//! \brief Return the total concentration of 'component'
//!
//! Only use this method if the mass is conserved for 'component'
//! \param component Index of the aqueous component (in the database)
scalar_t total_concentration_bc(index_t component) const {
specmicp_assert((component > 0 and aqueous_component_equation_type(component)
== AqueousComponentEquationType::MassConservation)
or
( water_equation_type() == WaterEquationType::MassConservation));
return m_fixed_values(component);
}
//! \brief Return the fixed saturation
//!
//! Only used if the saturation is fixed in the system
scalar_t fixed_saturation_bc() const {
specmicp_assert(water_equation_type() == WaterEquationType::FixedSaturation);
return m_fixed_values(0);
}
//! \brief Return the fixed activity of 'component'
//!
//! Only use this method if 'component' has a fixed activity constraint
//! \param component Index of the aqueous component (in the database)
scalar_t fixed_activity_bc(index_t component) const {
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedActivity);
return m_fixed_values(component);
}
//! \brief Return the fixed molality of 'component'
//!
//! Only use this method if 'component' has a fixed molality constraint
//! \param component Index of the aqueous component (in the database)
scalar_t fixed_molality_bc(index_t component) const {
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedMolality);
return m_fixed_values(component);
}
//! \brief Return the fixed fugacity value for 'component'
scalar_t fixed_fugacity_bc(index_t component) const {
specmicp_assert(aqueous_component_equation_type(component)
== AqueousComponentEquationType::FixedFugacity);
return m_fixed_values(component);
}
//! \brief Return the total concentration for the electron
scalar_t electron_total_concentration() const {
return 0.0;
}
// Gas
// ---
//! \brief Return the fugacity of 'gas'
//!
//! \param gas Index of the gas (in the database)
scalar_t gas_fugacity(index_t gas) const {
return m_second.gas_fugacity(gas);
}
//! \brief Return the concentration of 'gas' in the system
//!
//! \param gas Index of the gas (in the database)
scalar_t gas_concentration(index_t gas) const {
return m_second.gas_concentration(gas);
}
//! \brief Return true if gas is active
//!
//! \param gas Index of the gas (in the database)
bool is_active_gas(index_t gas) const {
return m_equations.active_gas[gas];
}
//! \brief Return the volume fraction (total saturation) of the gas phase
scalar_t volume_fraction_gas_phase() const noexcept {
return m_second.volume_fraction_gas;
}
// Sorbed species
// --------------
//! \brief Return the surface total concentration
const scalar_t& surface_total_concentration(index_t q) const {
return m_fixed_values(dof_surface(q));
}
//! \brief Return the sorption surface area
const scalar_t& sorption_surface_area() const {
return m_fixed_values(dof_surface_potential());
}
//! \brief Return true if 'sorbed' is an active species
//! \param sorbed Index of the sorbed species (in the database)
bool is_active_sorbed(index_t sorbed) const {
return m_equations.active_sorbed[sorbed];
}
//! \brief Return the molality of the sorbed species 'sorbed'
//! \param sorbed Index of the sorbed species (in the database)
const scalar_t& sorbed_species_concentration(index_t sorbed) const {
return m_second.sorbed_concentrations(sorbed);
}
// Pressure and temperature
// ------------------------
//! \brief Return the gas pressure
//!
//! This is a fixed value (1 atm)
scalar_t gas_total_pressure() const {
return 1.01325e5;
// Note : all pressure are in pascal, unit conversion is done for the concentration
}
//! \brief Return the temperature
//!
//! This is a fixed value (25°C)
scalar_t temperature() const noexcept {
return 273.15+25;
}
//! @}
// Solution
// =================
//! \brief Return the equilibrium state of the system, the Solution of the speciation problem
//!
//! Also checks that the secondary variables were correctly solved
//!
//! \param xtot the complete set of main variables
//! \param x the reduced set of main variables, only the variables with an equation
AdimensionalSystemSolution get_solution(Vector& xtot, const Vector& x);
//! \brief Return the equilibrium state of the system, the Solution of the speciation problem
//!
//! This version does not check anything.
//!
//! \param xtot the complete set of main variables
//! \param x the reduced set of main variables, only the variables with an equation
AdimensionalSystemSolution unsafe_get_solution(Vector& xtot, const Vector& x);
//! \name Secondary variables computation
//! \brief Algorithms to compute the secondary variables
//!
//! @{
//! \brief Compute the ionic strength
//!
//! \param x Vector of the main variables
void set_ionic_strength(const Vector& x);
//! \brief Compute the activity coefficients
//!
//! \param x Vector of the main variables
void compute_log_gamma(const Vector& x);
//! \brief Compute the secondary aqueous species molalities
//!
//! \param x Vector of the main variables
void set_secondary_concentration(const Vector& x);
//! \brief Compute the secondary variables
//!
//! \param x Vector of the main variables
void set_secondary_variables(const Vector& x);
//! \brief Compute the gas phase volume fraction
//!
//! \param x Vector of the main variables
void set_volume_fraction_gas_phase(const Vector& x);
//! \brief Compute the fugacity for all the gas
//!
//! \param x Vector of the main variables
void set_pressure_fugacity(const Vector& x);
//! \brief Compute the sorbed species concentrations
//!
//! \param x Vector of the main variables
void set_sorbed_concentrations(const Vector& x);
//! \brief This function is called at the beginning of each iteration by the solver
//!
//! It does the fixed-point iteration to compute the secondary variables
//! \param x Vector of the main variables
//! \param norm_residual Norm of the current residuals
bool hook_start_iteration(const Vector &x, scalar_t norm_residual);
//! @}
//! \brief A reasonable (well... maybe...) starting guess
//!
//! \param xtot Vector of the main variables
//! \sa #reasonable_starting_guess
void reasonable_starting_guess(Vector& xtot);
//! \brief A reasonable (maybe...) restarting guess
//!
//! \param xtot Vector of the main variables
//! \sa #reasonable_restarting_guess
void reasonable_restarting_guess(Vector& xtot);
//! \brief Set the units
void set_units(const units::UnitsSet& unit_set);
//! \brief Return true if the problem is a Box constrained VI
//! For BOX VI problem if needed
bool system_is_box_constrained() {
return m_equations.is_box_constrained;
}
//! \brief Flag an equation as a bo
bool is_box_vi(index_t ideq, scalar_t& upper_bound) {
upper_bound = m_equations.upper_bounds[ideq];
return m_equations.is_box_vi[ideq];
}
// Private Members and attributes
// ##############################
private:
void vector_debye_huckel_component(
scalar_t ionic_strength,
scalar_t sqrt_ionic
);
void vector_debye_huckel_aqueous(
scalar_t ionic_strength,
scalar_t sqrt_ionic
);
// scaling factors
void compute_scaling_factors();
//! \brief Sum of the aqueous molalities weighted by the stoichiometric coefficient of 'component'
scalar_t weigthed_sum_aqueous(index_t component) const;
//! \brief Sum of the sorbed concentration weighted by the stoichiometric coefficient of 'component'
scalar_t weigthed_sum_sorbed(index_t component) const;
//! \brief Sum of the sorbed concentration weighted by the stoichiometric coefficient of 'ssite'
scalar_t weigthed_sum_sorbed_ssite(index_t ssite) const;
//! \brief Sum of the sorbed concentration weighted by the charge of sorbed speces
scalar_t weigthed_sum_charge_sorbed() const;
//! \brief Sum of the mineral concentration weighted by the stoichiometric coefficient of 'component'
scalar_t weigthed_sum_mineral(const Vector& x, index_t component) const;
//! \brief Sum of the gas concentration weigthed by the stoichiometric coefficient of 'component'
scalar_t weigthed_sum_gas(index_t component) const;
//! \brief Derivative of the aqueous weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_aqueous
scalar_t diff_weigthed_sum_aqueous(index_t diff_component, index_t component) const;
//! \brief Derivative of the sorbed weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_sorbed
scalar_t diff_weigthed_sum_sorbed(index_t diff_component, index_t component) const;
//! \brief Derivative of the sorbed weigthed sum with respect to 'diff_ssite'
//! \sa weigthed_sum_sorbed
scalar_t diff_surface_weigthed_sum_sorbed(index_t diff_ssite, index_t component) const;
//! \brief Derivative of the sorbed weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_sorbed
scalar_t diff_weigthed_sum_sorbed_ssite(index_t diff_component, index_t ssite) const;
//! \brief Derivative of the sorbed weigthed sum with respect to 'diff_ssite'
//! \sa weigthed_sum_sorbed
scalar_t diff_surface_weigthed_sum_sorbed_ssite(index_t diff_ssite, index_t ssite) const;
//! \brief Return the square root in the second term of the surface potential in EDL model
scalar_t edl_sqrt_surface_potential() const;
//! \brief Derivative of the gas weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_sorbed
scalar_t diff_weigthed_sum_gas(index_t diff_component, index_t component) const;
// Jacobian
// ########
//! \brief Compute the water equation contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_water(Vector& x, Matrix& jacobian);
//! \brief Compute the aqueous components equation contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_aqueous_components(Vector& x, Matrix& jacobian);
//! \brief Compute the mineral equations contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_minerals(Vector& x, Matrix& jacobian);
//! \brief Compute the contribution of the surface sorption equation to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_surface(Vector& x, Matrix& jacobian);
//! \brief Compute the contribution of the surface potential equation to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_surface_potential(Vector& x, Matrix& jacobian);
//! \brief Compute the contribution of the electron equation to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void jacobian_electron(Vector& x, Matrix& jacobian);
// Equation numbering
// ##################
//! \brief Number the equations
//! \param constraints the constraints to apply to the system
//! \sa number_eq_aqueous_component
void number_eq(const AdimensionalSystemConstraints& constraints);
//! \brief Number the equations for the aqueous components
//! \param constraints the constraints to apply to the system
//! \param[in,out] neq the number of equations in the system
//! \sa number_eq
void number_eq_aqueous_component(
const AdimensionalSystemConstraints& constraints,
index_t& neq
);
// Setter
// ######
// main variables
// --------------
// they are set by the solver
// Secondary variables
// -------------------
//! \brief Return a reference to log_10(γ_i) where i is a component
scalar_t& log_gamma_component(index_t component) {
return m_second.loggamma(dof_component_gamma(component));
}
//! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species
scalar_t& log_gamma_secondary(index_t aqueous) {
return m_second.loggamma(dof_aqueous_gamma(aqueous));
}
//! \brief Return a reference to the ionic strength
scalar_t& ionic_strength() {
return m_second.ionic_strength;
}
// Attributes
// ##########
//! \struct SecondaryVariables
//! \brief Contains information about the secondary variables
struct SecondaryVariables {
scalar_t ionic_strength {0.0}; //!< The ionic Strength
scalar_t volume_fraction_gas {0.0}; //!< The gas saturation
scalar_t porosity {0.0}; //!< The porosity
Vector secondary_molalities; //!< The secondary molalities
Vector loggamma; //!< The log of activity coefficients
Vector gas_fugacity; //!< The gas fugacities
Vector gas_concentration; //!< The gas concentrations
Vector sorbed_concentrations; //!< The sorbed concentrations
//! \brief Initialization without a previous solution
SecondaryVariables(const RawDatabasePtr& data);
//! \brief Initialization with a previous solution
SecondaryVariables(const AdimensionalSystemSolution& previous_solution);
};
//! \struct IdEquations
//! \brief BookKeeper for the equations id and type
struct IdEquations {
bool use_water_pressure_model {false};
water_partial_pressure_f water_pressure_model {nullptr};
bool solve_pressure_model {false};
index_t nb_tot_variables {0};
index_t nb_free_variables {0};
index_t nb_complementarity_variables {0};
std::vector<index_t> ideq;
std::vector<index_t> component_equation_type;
std::vector<index_t> fixed_activity_species;
std::vector<index_t> nonactive_component {};
std::vector<bool> active_aqueous;
std::vector<bool> active_gas;
std::vector<bool> active_sorbed;
bool is_box_constrained {false};
std::vector<bool> is_box_vi {};
std::vector<scalar_t> upper_bounds {};
//! \brief Initialize the data structure
IdEquations(index_t nb_dofs, const RawDatabasePtr& data);
//! \brief Return the equation id
const index_t& id_equation(index_t dof) const {
return ideq[dof];
}
//! \brief Return a reference to the equation id
index_t& id_equation(index_t dof) {
return ideq[dof];
}
//! \brief Return a reference to the type of equation of 'dof'
index_t& type_equation(index_t dof) {
return component_equation_type[dof];
}
//! \brief Return the equation type for 'dof'
index_t type_equation(index_t dof) const {
return component_equation_type[dof];
}
//! \brief Return the related species for dof
//!
//! The exact meaning of this relation is dependant upon the equation type
index_t& related_species(index_t dof) {
return fixed_activity_species[dof];
}
//! \brief Add a non active component to the system
void add_non_active_component(index_t id) {
nonactive_component.push_back(id);
}
//! \brief Add an equation
void add_equation(index_t id, index_t *neq) {
ideq[id] = *neq;
++(*neq);
}
//! \brief Set the active flag of aqueous species 'id' to 'is_active'
void set_aqueous_active(index_t id, bool is_active) {
active_aqueous[id] = is_active;
}
//! \brief Set the active flag of gas 'id' to 'is_active'
void set_gas_active(index_t id, bool is_active) {
active_gas[id] = is_active;
}
//! \brief Set the active flag of sorbed species 'id' to 'is_active'
void set_sorbed_active(index_t id, bool is_active) {
active_sorbed[id] = is_active;
}
void activate_box_constrained(index_t nb_eq) {
is_box_constrained = true;
is_box_vi = std::vector<bool>(nb_eq, false);
upper_bounds = std::vector<scalar_t>(nb_eq, 0.0);
}
void set_box_constrained(index_t ideq_mineral, scalar_t upper_bound) {
is_box_vi[ideq_mineral] = true;
upper_bounds[ideq_mineral] = upper_bound;
}
};
bool not_in_linesearch {true};
scalar_t m_inert_volume_fraction;
scalar_t m_scaling_molar_volume {1.0};
scalar_t m_scaling_gas {1.0};
scalar_t m_scaling_molality {1.0};
Vector m_fixed_values {};
SecondaryVariables m_second;
IdEquations m_equations;
};
} // end namespace specmicp
#endif // SPECMICP_ADIMENSIONALSYSTEM_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_structs.hpp b/src/specmicp/adimensional/adimensional_system_structs.hpp
index 17c22b7..ae5dec1 100644
--- a/src/specmicp/adimensional/adimensional_system_structs.hpp
+++ b/src/specmicp/adimensional/adimensional_system_structs.hpp
@@ -1,441 +1,448 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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. *
============================================================================= */
#ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
#include "specmicp_common/types.hpp"
#include <vector>
//! \file adimensional_system_structs.hpp
//! \brief Options and constraints for the AdimensionalSystem
namespace specmicp {
//! \struct AdimensionalSystemOptions
//! \brief Options for the Adimensional Systems
//!
//! It is mainly about the secondary variables fixed-point iterations
struct SPECMICP_DLL_PUBLIC AdimensionalSystemOptions
{
//! Solve for non ideality
bool non_ideality;
//! Max iterations for the non ideality model
index_t non_ideality_max_iter;
//! Scaling the electron equation
scalar_t scaling_electron;
//! Tolerance for non ideality
scalar_t non_ideality_tolerance;
//! Under relaxation factor for the conservation of water
scalar_t under_relaxation_factor;
//! Log of the molality used to restart the computation
scalar_t restart_concentration;
//! Liquid water volume fraction to restart the computation
scalar_t restart_water_volume_fraction;
//! Log_10 of the molality for a new component
scalar_t new_component_concentration;
//! Factor to start the non-ideality computation
scalar_t start_non_ideality_computation;
//! Cutoff for including components in the computation
scalar_t cutoff_total_concentration;
AdimensionalSystemOptions(); // in adimensional_system_solver_structs.cpp
};
//! \enum AqueousComponentEquationType
//! \brief Type of an aqueous component equation
enum class AqueousComponentEquationType
{
NoEquation = no_equation, //!< Not an equation, component is not present in the system
MassConservation, //!< Mass balance
ChargeBalance, //!< M.B. replaced by charge balance
FixedFugacity, //!< M.B. replaced by a fixed fugacity equation
FixedActivity, //!< M.B. replaced by a fixed activity equation
FixedMolality //!< The molality of the species is fixed
};
//! \enum WaterEquationType
//! \brief The type of the equation solved for the water
enum class WaterEquationType
{
NoEquation = no_equation, //!< Amount of water is not solved
MassConservation, //!< Water is conserved
SaturatedSystem, //!< System is saturated
FixedSaturation //!< Saturation is fixed
};
//! \struct FixedFugacityConstraint
//! \brief Struct to contain information needed to solve a fix fugacity problem
struct FixedFugacityConstraint
{
index_t id_gas; //!< Index of the fixed-fugacity gas
index_t id_component; //!< Index of the corresponding component
scalar_t log_value; //!< Log_10 of the fugacity
//! \brief Constructor with initialization
FixedFugacityConstraint(index_t gas, index_t component, scalar_t logvalue):
id_gas(gas),
id_component(component),
log_value(logvalue)
{}
};
//! \struct FixedActivityConstraint
//! \brief Struct to contain information needed to solve a fix activity problem.
struct FixedActivityConstraint
{
index_t id_component; //!< Index of the fixed-activity component
scalar_t log_value; //!< Log_10 of the activity
//! \brief Constructor with initialization
FixedActivityConstraint(index_t component, scalar_t logvalue):
id_component(component),
log_value(logvalue)
{}
};
//! \struct FixedMolalityConstraint
//! \brief Contain information about a fixed molality constraint
struct FixedMolalityConstraint
{
index_t id_component; //!< Index of the component
scalar_t log_value; //!< Log_10 of the molality
//! \brief Constructor with initialization
FixedMolalityConstraint(index_t component, scalar_t logvalue):
id_component(component),
log_value(logvalue)
{}
};
//! \enum SurfaceEquationType
//! \brief The model for surface sorption
enum class SurfaceEquationType
{
NoEquation = no_equation, //!< Do not include surface sorption
Equilibrium, //!< Simple equilibrium model, no diffuse layer
EDL //!< Double layer diffuse model
};
//! \struct SurfaceConstraint
//! This struct contains the information to set-up the surface sorption model
struct SurfaceConstraint
{
SurfaceEquationType model_type; //!< The model to use
Vector total_ssite_concentrations; //!< The total concentration of sorption sites
scalar_t sorption_surface; //!< The sorption surface (area / volume material)
//! \brief By default, we don't include surface sorption in the computation
SurfaceConstraint() noexcept:
model_type(SurfaceEquationType::NoEquation),
sorption_surface(0.0)
{}
};
//! \brief Create a simple equilibrium surface model
inline SurfaceConstraint make_equilibrium_surface_model(
Vector surface_concentrations)
{
SurfaceConstraint model;
model.model_type = SurfaceEquationType::Equilibrium;
model.total_ssite_concentrations = surface_concentrations;
model.sorption_surface = 0.0;
return model;
}
//! \brief Create a EDL surface model
inline SurfaceConstraint make_EDL_surface_model(
Vector surface_concentrations,
scalar_t sorption_surface)
{
SurfaceConstraint model;
model.model_type = SurfaceEquationType::EDL;
model.total_ssite_concentrations = surface_concentrations;
model.sorption_surface = sorption_surface;
return model;
}
//! \enum ElectronEquationType the type of the equation for the electron
enum class ElectronEquationType
{
NoEquation = no_equation, //!< Do not compute the concentration equation of the electron
Equilibrium, //!< Set the concentration of electron to be 0
FixedpE //!< Activity of the electron is fixed
};
//! \enum MineralConstraintType
//! \brief Solid phase constraint type
enum class MineralConstraintType
{
Equilibrium, //!< Default - equilibrium phase
NoPrecipitation //!< No precipitation allowed
};
//! \brief Solid phase constraint
struct MineralConstraint
{
index_t id_mineral {no_species}; //!< Solid phase to be constrained
MineralConstraintType equation_type; //!< Type of the constraint
scalar_t param {0.0}; //!< Numerical value attached to the constraint
//! \brief Default constructor, no initialization
MineralConstraint() {}
//! \brief Constructor with initialization
MineralConstraint(
index_t mineral,
MineralConstraintType constraint,
scalar_t parameter
):
id_mineral(mineral),
equation_type(constraint),
param(parameter)
{}
};
//! \struct ElectronConstraint
//! \brief the constraint for the electron
struct ElectronConstraint
{
ElectronEquationType equation_type{ElectronEquationType::NoEquation}; //!< The equation type
scalar_t fixed_value; //!< The fixed value of pE if needed
index_t species {no_species}; //!< In case of fixed pE, this is the reaction to use
//! \brief By default we assume equilibrium
ElectronConstraint():
equation_type(ElectronEquationType::Equilibrium),
fixed_value(0.0)
{}
//! \brief When a value is provided, we assume that the pE is fixed
ElectronConstraint(scalar_t pe_value, scalar_t aqueous_species):
equation_type(ElectronEquationType::FixedpE),
fixed_value(pe_value),
species(aqueous_species)
{}
};
//! \brief Return the partial pressure as function of the liquid saturation
//!
//! This function must return the pressure in Pa
using water_partial_pressure_f = std::function<scalar_t (scalar_t)>;
//! \brief Constraints for the partial pressure model
struct WaterPartialPressureConstraint
{
bool use_partial_pressure_model; //!< True if we use partial pressure model
water_partial_pressure_f partial_pressure_model; //!< The model given by the user
//! \brief By default the partial pressure constraint is disabled
WaterPartialPressureConstraint():
use_partial_pressure_model(false),
partial_pressure_model(nullptr)
{}
//! \brief Enble the water partial pressure model
WaterPartialPressureConstraint(water_partial_pressure_f& model):
WaterPartialPressureConstraint()
{
if (model != nullptr)
{
use_partial_pressure_model = true;
partial_pressure_model = model;
}
}
};
//! \struct AdimensionalSystemConstraints
//! \brief Struct to contains the "Boundary conditions" for the AdimensionalSystem
//!
//! \ingroup specmicp_api
struct AdimensionalSystemConstraints
{
//! Total concentrations
Vector total_concentrations;
// Water
//! Water equation
WaterEquationType water_equation {WaterEquationType::MassConservation};
//! The partial pressure of water (if any)
WaterPartialPressureConstraint water_partial_pressure {};
//! \brief A parameter for the water equation
//!
//! It is the saturation if the system has a fixed saturation,
//! or invalid otherwise
scalar_t water_parameter {-1.0};
//! The equation for this component is replace by the charge balance
index_t charge_keeper{no_species};
//! Contains information about fixed fugacity constraints
std::vector<FixedFugacityConstraint> fixed_fugacity_cs {};
//! Contains information about fixed activity constraints
std::vector<FixedActivityConstraint> fixed_activity_cs {};
//! Contains information about fixed molality constraints
std::vector<FixedMolalityConstraint> fixed_molality_cs {};
//! Additional constraints for solid phases
std::vector<MineralConstraint> mineral_constraints {};
//! \brief Volume fraction of inert solid
//!
//! (inert in the equilibrium computation)
scalar_t inert_volume_fraction{ 0.0 };
SurfaceConstraint surface_model{}; //!< Surface sorption model
ElectronConstraint electron_constraint{}; //!< constraint for the electron
+ bool immobile_disabled {false}; //!< True if immobile species are not solved
+
//! \brief Default constructor
AdimensionalSystemConstraints():
total_concentrations()
{}
//! \brief Constructor, initialize total concentrations
AdimensionalSystemConstraints(const Vector& total_concs):
total_concentrations(total_concs)
{}
//! \brief Set the total concentrations
void set_total_concentrations(const Vector& total_concs) {total_concentrations = total_concs;}
//! \brief Enable the conservation of water
void enable_conservation_water() noexcept {water_equation = WaterEquationType::MassConservation;}
//! \brief Disable the conservation of water
void disable_conservation_water() noexcept {water_equation = WaterEquationType::NoEquation;}
//! \brief The system is saturated
void set_saturated_system() noexcept {water_equation = WaterEquationType::SaturatedSystem;}
//! \brief Fixed the saturation in the system
void set_fixed_saturation(scalar_t saturation) noexcept {
water_equation = WaterEquationType::FixedSaturation;
water_parameter = saturation;
}
//! \brief Disable the surface sorption model
void disable_surface_model() noexcept {surface_model = SurfaceConstraint();}
//! \brief Create a simple equilibrium surface model
void set_equilibrium_surface_model(Vector ssites_concentrations) {
surface_model = make_equilibrium_surface_model(ssites_concentrations);
}
//! \brief Create an EDL surface model
void set_EDL_surface_model(Vector ssites_concentrations, scalar_t sorption_surface) {
surface_model = make_EDL_surface_model(ssites_concentrations, sorption_surface);
}
//! \brief Set the charge keeper to 'component'
//!
//! \param component Index of the component (in the database)
void set_charge_keeper(index_t component) noexcept {
charge_keeper = component;
}
//! \brief Add a fixed fugacity gas condition
//!
//! \param constraint struct containing the information about a fixed-fugacity constraint
void add_fixed_fugacity_gas(const FixedFugacityConstraint& constraint) {
fixed_fugacity_cs.push_back(constraint);
}
//! \brief Add a fixed fugacity gas condition
//!
//! \param gas Index of the gas (in the database)
//! \param component Index of the corresponding component (in the database)
//! \param logvalue Log_10 of the fugacity
void add_fixed_fugacity_gas(index_t gas, index_t component, scalar_t logvalue) {
fixed_fugacity_cs.emplace_back(FixedFugacityConstraint(gas, component, logvalue));
}
//! \brief Add a fixed activity component condition
//!
//! \param constraint struct containing the information about a fixed-activity constraint
void add_fixed_activity_component(const FixedActivityConstraint& constraint) {
fixed_activity_cs.push_back(constraint);
}
//! \brief Add a fixed activity condition for a component
//!
//! \param component Index of the corresponding component (in the database)
//! \param log_value Log_10 of the activity
void add_fixed_activity_component(index_t component, scalar_t log_value) {
fixed_activity_cs.emplace_back(FixedActivityConstraint(component, log_value));
}
//! \brief Add a fixed molality condition for a component
//!
//! \param component Index of the corresponding component (in the database)
//! \param log_value Log_10 of the molality
void add_fixed_molality_component(index_t component, scalar_t log_value) {
fixed_molality_cs.emplace_back(FixedMolalityConstraint(component, log_value));
}
//! \brief Set the inert volume fraction
//!
//! The volume fraction of the inert phase is used to offset the saturation.
//! This inert phase may correspond to aggregates or solid phases governed by kinetics.
//!
//! \param value volume fraction of the inert phase
void set_inert_volume_fraction(scalar_t value) noexcept {
inert_volume_fraction = value;
}
//! \brief Set a partial pressure model for water
//!
//! The partial pressure model for water is a function taking the saturation
//! and returning the partial pressure of water
void set_water_partial_pressure_model(water_partial_pressure_f model)
{
water_partial_pressure.use_partial_pressure_model = true;
water_partial_pressure.partial_pressure_model = model;
}
//! \brief Set an upper bound for the volume fraction of the mineral
//!
//! Emulate a no precipitation condition
void set_mineral_upper_bound(index_t mineral, scalar_t upper_volume_fraction)
{
mineral_constraints.emplace_back(mineral,
MineralConstraintType::NoPrecipitation,
upper_volume_fraction);
}
+ //! \brief Disable the immobile species (solid phases and sorbed species)
+ void disable_immobile_species() noexcept {
+ immobile_disabled = true;
+ }
+
// //! \brief Set the system at a fixed pE
// void set_fixed_pe(scalar_t pe_value, index_t aqueous_species) noexcept {
// electron_constraint = ElectronConstraint(pe_value, aqueous_species);
// }
};
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
diff --git a/src/specmicp/io/configuration.cpp b/src/specmicp/io/configuration.cpp
index 374fea4..1e8495e 100644
--- a/src/specmicp/io/configuration.cpp
+++ b/src/specmicp/io/configuration.cpp
@@ -1,842 +1,856 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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 "configuration.hpp"
#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
#include "specmicp/adimensional/config_default_options_solver.h"
#include "specmicp_common/io/safe_config.hpp"
#include "specmicp_common/physics/units.hpp"
#include "specmicp_common/physics/io/configuration.hpp"
#include "specmicp_common/physics/laws.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "specmicp_database/database.hpp"
#include "specmicp_common/io/config_yaml_sections.h"
#include <unordered_map>
namespace specmicp {
namespace io {
// Additional declarations
// =======================
//! \brief Configure a fixed fugacity component
//!
//! \deprecated
void configure_specmicp_constraints_fixed_activity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_activity,
const RawDatabasePtr& raw_db
);
//! \brief Configure a fixed fugacity component
//!
//! \deprecated
void configure_specmicp_constraints_fixed_fugacity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_fugacity,
const RawDatabasePtr& raw_db
);
//! \brief Configure a fixed molality constraint
//!
//! \deprecated
void configure_specmicp_constraints_fixed_molality(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_molality,
const RawDatabasePtr& raw_db
);
//! \brief Configure a fixed fugacity component
void configure_specmicp_constraints_fixed_activity(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_activity
);
//! \brief Configure a fixed fugacity component
void configure_specmicp_constraints_fixed_fugacity(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_fugacity
);
//! \brief Configure a fixed molality constraint
void configure_specmicp_constraints_fixed_molality(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_molality
);
// Implementation
// ==============
void configure_specmicp_options(
AdimensionalSystemSolverOptions& options,
const units::UnitsSet& the_units,
YAMLConfigHandle&& conf
)
{
micpsolver::MiCPSolverOptions& solver = options.solver_options;
conf.set_if_attribute_exists<scalar_t>(
solver.fvectol, SPC_CF_S_SPECMICP_A_RES_TOL, 0.0);
conf.set_if_attribute_exists<scalar_t>(
solver.steptol, SPC_CF_S_SPECMICP_A_STEP_TOL, 0.0);
conf.set_if_attribute_exists<index_t>(
solver.max_iter, SPC_CF_S_SPECMICP_A_MAX_ITER, 0);
conf.set_if_attribute_exists<scalar_t>(
solver.maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH, 0.0);
conf.set_if_attribute_exists<index_t>(
solver.maxiter_maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER, 0);
conf.set_if_attribute_exists<bool>(
solver.use_scaling, SPC_CF_S_SPECMICP_A_SCALING);
conf.set_if_attribute_exists<bool>(
solver.non_monotone_linesearch, SPC_CF_S_SPECMICP_A_NONMONOTONE);
conf.set_if_attribute_exists<scalar_t>(
solver.factor_descent_condition, SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION);
conf.set_if_attribute_exists<scalar_t>(
solver.condition_limit, SPC_CF_S_SPECMICP_A_COND_CHECK);
conf.set_if_attribute_exists<scalar_t>(
solver.threshold_cycling_linesearch, SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH);
AdimensionalSystemOptions& system = options.system_options;
conf.set_if_attribute_exists<scalar_t>(
system.non_ideality_tolerance, SPC_CF_S_SPECMICP_A_NONIDEAL_TOL, 0.0);
conf.set_if_attribute_exists<index_t>(
system.non_ideality_max_iter, SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER, 0);
conf.set_if_attribute_exists<scalar_t>(
system.cutoff_total_concentration, SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC, 0.0);
conf.set_if_attribute_exists<scalar_t>(
system.restart_concentration, SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION);
conf.set_if_attribute_exists<scalar_t>(
system.restart_water_volume_fraction, SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC, 0.0);
conf.set_if_attribute_exists<scalar_t>(
system.under_relaxation_factor, SPC_CF_S_SPECMICP_A_UNDER_RELAXATION, 0.0, 1.0);
options.units_set = the_units;
}
// ReactantBox
// ============
namespace internal {
// helper functions
static void configure_formulation(
ReactantBox& reactant_box,
YAMLConfigHandle&& configuration
);
static void configure_formulation_aqueous(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_aqueous
);
static void configure_formulation_solid(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_aqueous
);
static void configure_constraints(
ReactantBox& reactant_box,
YAMLConfigHandle&& configuration
);
static void configure_constraints_fixed_activity(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_act
);
static void configure_constraints_fixed_fugacity(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_act
);
static void configure_constraints_fixed_molality(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_mol
);
} // end namespace internal
void SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box(
ReactantBox& reactant_box,
YAMLConfigHandle&& configuration
)
{
internal::configure_formulation(
reactant_box,
configuration.get_section(SPC_CF_S_FORMULATION)
);
internal::configure_constraints(
reactant_box,
configuration.get_section(SPC_CF_S_CONSTRAINTS)
);
}
ReactantBox SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box(
RawDatabasePtr raw_db,
const units::UnitsSet& the_units,
YAMLConfigHandle&& configuration
)
{
ReactantBox reactant_box(raw_db, the_units);
configure_specmicp_reactant_box(reactant_box, std::move(configuration));
return reactant_box;
}
namespace internal {
// helper functions
static void configure_formulation(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_formulation)
{
// solution
auto cf_solution = cf_formulation.get_section(SPC_CF_S_FORMULATION_A_SOLUTION);
const auto amount =
cf_solution.get_required_attribute<double>(SPC_CF_S_FORMULATION_A_AMOUNT);
const auto unit =
cf_solution.get_required_attribute<std::string>(SPC_CF_S_FORMULATION_A_UNIT);
reactant_box.set_solution(amount, unit);
if (cf_formulation.has_section(SPC_CF_S_FORMULATION_A_AQUEOUS))
{
configure_formulation_aqueous(
reactant_box,
cf_formulation.get_section(SPC_CF_S_FORMULATION_A_AQUEOUS)
);
}
if (cf_formulation.has_section(SPC_CF_S_FORMULATION_A_MINERALS))
{
configure_formulation_solid(
reactant_box,
cf_formulation.get_section(SPC_CF_S_FORMULATION_A_MINERALS)
);
}
}
static void configure_formulation_aqueous(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_aqueous
)
{
if (not cf_aqueous.is_sequence()) {
cf_aqueous.report_error(
YAMLConfigError::ListExpected,
"Aqueous species must be provide as a list of triplet"
" {label amount, unit}"
);
}
for (auto ind: RangeIterator<uindex_t>(cf_aqueous.size()))
{
auto cf = cf_aqueous.get_section(ind);
const auto label =
cf.get_required_attribute<std::string>(SPC_CF_S_FORMULATION_A_LABEL);
const auto value =
cf.get_required_attribute<double>(SPC_CF_S_FORMULATION_A_AMOUNT);
const auto unit =
cf.get_required_attribute<std::string>(SPC_CF_S_FORMULATION_A_UNIT);
reactant_box.add_aqueous_species(label, value, unit);
}
}
static void configure_formulation_solid(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_solid
)
{
if (not cf_solid.is_sequence()) {
cf_solid.report_error(
YAMLConfigError::ListExpected,
"Solid species must be provide as a list of triplet"
" {label amount, unit}"
);
}
for (auto ind: RangeIterator<uindex_t>(cf_solid.size()))
{
auto cf = cf_solid.get_section(ind);
const auto label = cf.get_required_attribute<
std::string>(SPC_CF_S_FORMULATION_A_LABEL);
const auto value = cf.get_required_attribute<
double>(SPC_CF_S_FORMULATION_A_AMOUNT);
const auto unit = cf.get_required_attribute<
std::string>(SPC_CF_S_FORMULATION_A_UNIT);
reactant_box.add_solid_phase(label, value, unit);
}
}
static void configure_constraints(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf)
{
// find the constraint to apply
bool has_saturated = false;
bool has_fixed_saturation = false;
bool has_water_conservation = false;
bool does_conservation_water = true;
scalar_t water_param = -1.0;
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM))
{
has_saturated = cf.get_attribute<bool>(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM);
}
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_FIXEDSATURATION))
{
water_param = cf.get_attribute<scalar_t>(SPC_CF_S_CONSTRAINTS_A_FIXEDSATURATION, 0.0, 1.0);
if (water_param == 1.0)
{
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) {
cf.report_error(YAMLConfigError::InvalidArgument,
"SaturatedSystem and a Fixed saturation of 1.0 are both set."
"Provide only one of the options");
}
has_saturated = true;
}
else
has_fixed_saturation = true;
}
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER))
{
has_water_conservation = true;
does_conservation_water = cf.get_attribute<bool>(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER);
}
else
{
has_water_conservation = false;
does_conservation_water = true;
}
if ( (int) has_saturated + (int) has_fixed_saturation + (int) has_water_conservation > 1)
{
cf.report_error(YAMLConfigError::InvalidArgument,
"Incompatible options, choose only one of : "
"Saturated system, fixed saturation or conservation water");
}
if (has_saturated) {
reactant_box.set_saturated_system();
}
else if (has_fixed_saturation) {
reactant_box.set_fixed_saturation(water_param);
}
else if (not does_conservation_water) {
reactant_box.disable_conservation_water();
}
// Charge keeper
// -------------
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER))
{
auto label = cf.get_attribute<
std::string>(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER);
reactant_box.set_charge_keeper(label);
}
// Fixed activity
// ---------------
if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY))
{
configure_constraints_fixed_activity(
reactant_box,
cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY)
);
}
// Fixed fugacity
// --------------
if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY))
{
configure_constraints_fixed_fugacity(
reactant_box,
cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY)
);
}
// Fixed molality
// ---------------
if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY))
{
configure_constraints_fixed_molality(
reactant_box,
cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY)
);
}
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION))
{
reactant_box.set_inert_volume_fraction(
cf.get_attribute<double>(SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION)
);
}
+
+ // Immobile species
+ // ----------------
+ if (cf.get_optional_attribute<bool>(SPC_CF_S_CONSTRAINTS_A_DISABLEIMMOBILE, false))
+ {
+ reactant_box.disable_immobile_species();
+ }
}
static void configure_constraints_fixed_activity(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_act
)
{
if (not cf_f_act.is_sequence()) {
cf_f_act.report_error(
YAMLConfigError::ListExpected,
"Fixed activity component must be provide as a list of"
" triplet {label amount, unit}"
);
}
for (auto ind: RangeIterator<uindex_t>(cf_f_act.size()))
{
auto cf = cf_f_act.get_section(ind);
const auto label = cf.get_required_attribute<
std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL);
double value {0.0}; // value can be provided as value or log_value
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT))
{
value = cf.get_attribute<double>(SPC_CF_S_CONSTRAINTS_A_AMOUNT);
}
else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG))
{
const auto log_val = cf.get_attribute<
double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG);
value = std::pow(10.0, log_val);
}
else
{
cf.report_error(
YAMLConfigError::MissingRequiredAttribute,
"Either amount or log amount must be provided."
);
}
reactant_box.add_fixed_activity_component(label, value);
}
}
static void configure_constraints_fixed_fugacity(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_fug
)
{
if (not cf_f_fug.is_sequence()) {
cf_f_fug.report_error(
YAMLConfigError::ListExpected,
"Fixed fugacity gas must be provide as a list of"
" triplet {label amount, unit}"
);
}
for (auto ind: RangeIterator<uindex_t>(cf_f_fug.size()))
{
auto cf = cf_f_fug.get_section(ind);
// both component and gas must be provided
const auto label_c = cf.get_required_attribute<
std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT);
const auto label_g = cf.get_required_attribute<
std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL_GAS);
double value {0.0}; // value can be provided as value or log_value
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT))
{
value = cf.get_attribute<double>(SPC_CF_S_CONSTRAINTS_A_AMOUNT);
}
else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG))
{
const auto log_val = cf.get_attribute<
double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG);
value = std::pow(10.0, log_val);
}
else
{
cf.report_error(
YAMLConfigError::MissingRequiredAttribute,
"Either amount or log amount must be provided."
);
}
reactant_box.add_fixed_fugacity_gas(label_g, label_c, value);
}
}
static void configure_constraints_fixed_molality(
ReactantBox& reactant_box,
YAMLConfigHandle&& cf_f_mol
)
{
if (not cf_f_mol.is_sequence()) {
cf_f_mol.report_error(
YAMLConfigError::ListExpected,
"Fixed molality component must be provide as a list of"
" triplet {label amount, unit}"
);
}
for (auto ind: RangeIterator<uindex_t>(cf_f_mol.size()))
{
auto cf = cf_f_mol.get_section(ind);
const auto label = cf.get_required_attribute<
std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL);
double value {0.0}; // value can be provided as value or log_value
if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT))
{
value = cf.get_attribute<double>(SPC_CF_S_CONSTRAINTS_A_AMOUNT);
}
else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG))
{
const auto log_val = cf.get_attribute<
double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG);
value = std::pow(10.0, log_val);
}
else
{
cf.report_error(
YAMLConfigError::MissingRequiredAttribute,
"Either amount or log amount must be provided."
);
}
reactant_box.add_fixed_molality_component(label, value);
}
}
} // end namespace internal
// Constraints
// ===========
// new interface
// -------------
void configure_specmicp_constraints(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints
)
{
// equation for water
// ==================
// enable conservation of water
if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER))
{
auto tmp_water = conf_constraints.get_attribute<bool>(
SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER);
if (tmp_water)
{
// only one constraints
if (conf_constraints.get_optional_attribute<bool>(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM, false))
{
conf_constraints.report_error(YAMLConfigError::InvalidArgument,
"Attributes " SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER
"and " SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM
" cannot be set to true at the same time");
}
constraints.enable_conservation_water();
}
else
constraints.disable_conservation_water();
}
// or enable saturated system
if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM))
{
if (conf_constraints.get_attribute<bool>(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM))
{
constraints.set_saturated_system();
}
}
conf_constraints.set_if_attribute_exists(constraints.inert_volume_fraction,
SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION
);
// Aqueous components
// ==================
// charge keeper
// -------------
if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER))
{
const auto label = conf_constraints.get_attribute<std::string>(
SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER);
const index_t id = raw_db->get_id_component(label);
if (id == no_species)
{
conf_constraints.report_error(
YAMLConfigError::InvalidArgument,
"Species " + label + " does not exist in the database."
" It can't be the charge keeper."
);
}
constraints.set_charge_keeper(id);
}
// Fixed activity
// ---------------
if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY))
{
configure_specmicp_constraints_fixed_activity(
constraints,
raw_db,
conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY)
);
}
// Fixed fugacity
// --------------
if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY))
{
configure_specmicp_constraints_fixed_fugacity(
constraints,
raw_db,
conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY)
);
}
// Fixed activity
// ---------------
if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY))
{
configure_specmicp_constraints_fixed_molality(
constraints,
raw_db,
conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY)
);
}
// Surface model
// =============
if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL))
{
const bool is_enabled = conf_constraints.get_attribute<bool>(
SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL
);
if (is_enabled)
{
const auto value = conf_constraints.get_required_attribute<scalar_t>(
SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION);
Vector x;
x.resize(1);
x(0) = value;
constraints.set_equilibrium_surface_model(x);
}
else
{
constraints.disable_surface_model();
}
}
+ // Immobile concentration
+ // ======================
+ if (conf_constraints.get_optional_attribute<bool>(
+ SPC_CF_S_CONSTRAINTS_A_DISABLEIMMOBILE, false))
+ {
+ constraints.disable_immobile_species();
+ }
}
// helpers functions
//! \brief Return the id of a component
index_t get_id_component(
const database::DataContainer* const raw_db,
YAMLConfigHandle& conf,
const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL
)
{
const auto label = conf.get_required_attribute<std::string>(attr_key);
const index_t id = raw_db->get_id_component(label);
if (id == no_species)
{
conf.report_error(
YAMLConfigError::InvalidArgument,
"'" + label + "' is not a valid component."
);
}
return id;
}
//! \brief Return the id of a gas
index_t get_id_gas(
const database::DataContainer* const raw_db,
YAMLConfigHandle& conf,
const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL_GAS
)
{
const auto label = conf.get_required_attribute<std::string>(attr_key);
const index_t id = raw_db->get_id_gas(label);
if (id == no_species)
{
conf.report_error(
YAMLConfigError::InvalidArgument,
"'" + label + "' is not a valid gas."
);
}
return id;
}
//! \brief Return a value in log10 for fixed fugacity
scalar_t get_value(YAMLConfigHandle& conf)
{
scalar_t value = 0;
if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) {
value = std::log10(conf.get_attribute<scalar_t>(SPC_CF_S_CONSTRAINTS_A_AMOUNT));
}
else if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) {
value = conf.get_attribute<scalar_t>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG);
}
else {
conf.report_error(
YAMLConfigError::InvalidArgument,
"Expected one of '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '"
SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "'."
);
}
return value;
}
void configure_specmicp_constraints_fixed_activity(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_activity
)
{
for (uindex_t ind=0; ind<conf_constraints_fixed_activity.size(); ++ind)
{
YAMLConfigHandle conf = conf_constraints_fixed_activity.get_section(ind);
const index_t id = get_id_component(raw_db, conf);
const scalar_t value = get_value(conf);
constraints.add_fixed_activity_component(id, value);
}
}
void configure_specmicp_constraints_fixed_fugacity(
AdimensionalSystemConstraints& constraints,
const database::DataContainer* const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_fugacity
)
{
for (uindex_t ind=0; ind<conf_constraints_fixed_fugacity.size(); ++ind)
{
YAMLConfigHandle conf = conf_constraints_fixed_fugacity.get_section(ind);
const index_t id_c = get_id_component(raw_db, conf, SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT);
const index_t id_g = get_id_gas(raw_db, conf, SPC_CF_S_CONSTRAINTS_A_LABEL_GAS);
const scalar_t value = get_value(conf);
constraints.add_fixed_fugacity_gas(id_g, id_c, value);
}
}
void configure_specmicp_constraints_fixed_molality(
AdimensionalSystemConstraints& constraints,
const database::DataContainer * const raw_db,
YAMLConfigHandle&& conf_constraints_fixed_molality
)
{
for (uindex_t ind=0; ind<conf_constraints_fixed_molality.size(); ++ind)
{
YAMLConfigHandle conf = conf_constraints_fixed_molality.get_section(ind);
const index_t id = get_id_component(raw_db, conf);
const scalar_t value = get_value(conf);
constraints.add_fixed_molality_component(id, value);
}
}
void configure_smart_solver_initialization(
SmartAdimSolver& solver,
YAMLConfigHandle&& conf_init
)
{
if (conf_init.has_attribute(SPC_CF_S_INITIALIZATION_A_SOLUTION))
{
solver.set_init_volfrac_water(
conf_init.get_attribute<scalar_t>(SPC_CF_S_INITIALIZATION_A_SOLUTION));
}
if (conf_init.has_section(SPC_CF_S_INITIALIZATION_A_AQUEOUS))
{
auto scf = conf_init.get_section(SPC_CF_S_INITIALIZATION_A_AQUEOUS);
if (not scf.is_map())
{
scf.report_error(YAMLConfigError::MapExpected,
"A map {component, value} is expected.");
}
if (scf.has_attribute("all"))
{
solver.set_init_molality(scf.get_attribute<scalar_t>("all"));
}
std::unordered_map<std::string, scalar_t> tmap;
tmap.reserve(scf.size());
for (auto it=scf.map_begin(); it!=scf.map_end(); ++it)
{
auto name = (*it).first;
if (name != "all")
{
tmap.insert({name, scf.get_attribute<scalar_t>(name)});
}
}
if (not tmap.empty())
{
solver.set_init_molality(tmap);
}
}
if (conf_init.has_section(SPC_CF_S_INITIALIZATION_A_MINERALS))
{
auto scf = conf_init.get_section(SPC_CF_S_INITIALIZATION_A_MINERALS);
if (not scf.is_map())
{
scf.report_error(YAMLConfigError::MapExpected,
"A map {component, value} is expected.");
}
std::unordered_map<std::string, scalar_t> tmap;
tmap.reserve(scf.size());
for (auto it=scf.map_begin(); it!=scf.map_end(); ++it)
{
auto name = (*it).first;
tmap.insert({name, scf.get_attribute<scalar_t>(name)});
}
if (not tmap.empty())
{
solver.set_init_volfrac_mineral(tmap);
}
}
}
} //end namespace io
} //end namespace specmicp
diff --git a/src/specmicp/io/configuration.hpp b/src/specmicp/io/configuration.hpp
index 66a669c..ef7ff25 100644
--- a/src/specmicp/io/configuration.hpp
+++ b/src/specmicp/io/configuration.hpp
@@ -1,115 +1,115 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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. *
============================================================================= */
#ifndef SPECMICP_SPECMICP_IO_HPP
#define SPECMICP_SPECMICP_IO_HPP
#include "specmicp_common/types.hpp"
#include <memory>
//! \file specmicp/io/configuration.hpp
//! \brief configuration for specmicp adim system
namespace specmicp {
struct AdimensionalSystemSolverOptions;
struct Formulation;
struct AdimensionalSystemConstraints;
class ReactantBox;
class SmartAdimSolver;
namespace database {
struct DataContainer;
}
using RawDatabasePtr = std::shared_ptr<database::DataContainer>;
namespace units {
struct UnitsSet;
} //end namespace units
namespace io {
class YAMLConfigHandle;
/*! \brief Configure specmicp solver options
Conserve previous values if not given.
Configuration must be of the form :
\code
opts1: val1
opts2: val2
....
\endcode
\warning units must be set elsewhere !
*/
void SPECMICP_DLL_PUBLIC configure_specmicp_options(
AdimensionalSystemSolverOptions& options,
const units::UnitsSet& the_units,
YAMLConfigHandle&& configuration
);
//! \brief Configure a reactant box from the configuration
ReactantBox SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box(
RawDatabasePtr raw_db,
const units::UnitsSet& the_units,
YAMLConfigHandle&& configuration
);
//! \brief Configure a reactant box from the configuration
void SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box(
ReactantBox& reactant_box,
YAMLConfigHandle&& configuration
);
//! \brief Configure specmicp constraints
void SPECMICP_DLL_PUBLIC configure_specmicp_constraints(
AdimensionalSystemConstraints& constraints,
const database::DataContainer* const raw_db,
YAMLConfigHandle&& conf_constraints
);
//! \brief Configure AdimSmartSolver initialization
void SPECMICP_DLL_PUBLIC configure_smart_solver_initialization(
SmartAdimSolver& solver,
YAMLConfigHandle&& conf_init
);
} //end namespace io
} //end namespace specmicp
#endif // SPECMICP_SPECMICP_IO_HPP
diff --git a/src/specmicp/problem_solver/reactant_box.cpp b/src/specmicp/problem_solver/reactant_box.cpp
index dd105c3..71f337f 100644
--- a/src/specmicp/problem_solver/reactant_box.cpp
+++ b/src/specmicp/problem_solver/reactant_box.cpp
@@ -1,715 +1,722 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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 {
+ bool immobile_disabled {false};
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);
}
+void ReactantBox::disable_immobile_species()
+{
+ m_impl->immobile_disabled = true;
+}
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);
-
+ if ( m_impl->immobile_disabled ) {
+ constraints.disable_immobile_species();
+ }
// 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
diff --git a/src/specmicp/problem_solver/reactant_box.hpp b/src/specmicp/problem_solver/reactant_box.hpp
index 27c58f5..7e2e800 100644
--- a/src/specmicp/problem_solver/reactant_box.hpp
+++ b/src/specmicp/problem_solver/reactant_box.hpp
@@ -1,204 +1,206 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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. *
============================================================================= */
#ifndef SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
#define SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
//! \file problem_solver/reactant_box.hpp
//! \brief Handle initial composition, with units
#include "specmicp_common/types.hpp"
#include "specmicp_common/pimpl_ptr.hpp"
#include "specmicp_common/physics/units.hpp"
#include "specmicp_database/database_holder.hpp"
namespace specmicp {
struct AdimensionalSystemConstraints;
using water_partial_pressure_f = std::function<scalar_t (scalar_t)>;
//! \brief A set of reactants
//!
//! Take care of the units
class SPECMICP_DLL_PUBLIC ReactantBox :
public units::UnitBaseClass,
public database::DatabaseHolder
{
public:
//! \brief Constructor
ReactantBox(
std::shared_ptr<database::DataContainer> raw_data,
const units::UnitsSet& units_set
);
//! \brief Destructor
~ReactantBox();
//! \brief Move constructor
ReactantBox(ReactantBox&& other);
//! \brief Move constructor
ReactantBox& operator=(ReactantBox&& other);
// Solution
// =========
//! \brief Set the solution amount
void set_solution(scalar_t value, std::string str_unit);
//! \brief Set the solution amount
void set_solution(
scalar_t value,
const units::AmountUnit& unit
);
// Aqueous
// =======
//! \brief Add an aqueous species
void add_aqueous_species(
std::string name,
scalar_t value,
std::string str_unit
);
//! \brief Add an aqueous species
void add_aqueous_species(
std::string name,
scalar_t value,
const units::AmountUnit& unit
);
//! \brief Set the amount of an aqueous species
void set_aqueous_species(
std::string name,
scalar_t value,
std::string str_unit
);
//! \brief Set the amount an aqueous species
void set_aqueous_species(
std::string name,
scalar_t value,
const units::AmountUnit& unit
);
// Solid
// =====
//! \brief Add a solid phase
void add_solid_phase(
std::string name,
scalar_t value,
std::string str_unit
);
//! \brief Add a solid phase
void add_solid_phase(
std::string name,
scalar_t value,
const units::AmountUnit& unit
);
//! \brief Set the amount of a solid phase
void set_solid_phase(
std::string name,
scalar_t value,
std::string str_unit
);
//! \brief Set the amount of a solid phase
void set_solid_phase(
std::string name,
scalar_t value,
const units::AmountUnit& unit
);
// Total concentration
// ===================
//! \brief add a component without corresponding species
//!
//! It is useful if the database must be modified
void add_component(const std::string& name);
//! \brief Return the total concentration vector
Vector get_total_concentration(bool modify_db=false) const;
// Constraints
// ===========
//! \brief Set the charge keeper
void set_charge_keeper(std::string charge_keeper);
//! \brief Set the inert volume fraction (0 by default)
void set_inert_volume_fraction(scalar_t inert_volume_fraction);
//! \brief The system is saturated
void set_saturated_system();
//! \brief The saturation of the system is fixed
void set_fixed_saturation(scalar_t saturation);
//! \brief Disable the water conservation equation
//!
//! \warning probably not a good idead
void disable_conservation_water();
//! \brief Set the water partial pressure model
void set_water_partial_pressure_model(water_partial_pressure_f h2o_pressure_model);
//! \brief Add a fixed fugacity gas
void add_fixed_fugacity_gas(
std::string gas,
std::string component,
scalar_t fugacity
);
//! \brief Add a fixed activity component
void add_fixed_activity_component(
std::string component,
scalar_t activity
);
//! \brief Add a fixed molality component
//!
//! Molality is in mol/kg !
void add_fixed_molality_component(
std::string component,
scalar_t molality
);
-
//! \brief Set an upper bound for a mineral
//!
//! It's an upper bound on the volume fraction
void set_mineral_upper_bound(
std::string mineral,
scalar_t max_volume_fraction
);
+ //! \brief Disable the immobile species while solving the system
+ void disable_immobile_species();
+
//! \brief Return a set of constraints to be passed to the specmicp solver
AdimensionalSystemConstraints get_constraints(bool modify_db=true) const;
private:
struct SPECMICP_DLL_LOCAL ReactantBoxImpl;
utils::pimpl_ptr<ReactantBoxImpl> m_impl;
};
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
diff --git a/src/specmicp_common/io/config_yaml_sections.h b/src/specmicp_common/io/config_yaml_sections.h
index 711b3d9..8cebd94 100644
--- a/src/specmicp_common/io/config_yaml_sections.h
+++ b/src/specmicp_common/io/config_yaml_sections.h
@@ -1,826 +1,827 @@
/* =============================================================================
- Copyright (c) 2014 - 2016
- F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
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. *
============================================================================= */
#ifndef SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H
#define SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H
//! \file io/config_yaml_sections.h
//! \brief Define sections in the YAML config files
//!
//! Name are built as follow
//! PREFIX_[X_NAME]n
//!
//! where PREFIX is SPC_CF (specmicp conf)
//! X is either A (for attribute), v (for value), S (for section) or SS (for subsection)
//! and NAME is the name used in the conf file
//! \def SPC_CF_NAME_PATH_ROOT
//! \brief Name of the first level
#define SPC_CF_NAME_PATH_ROOT "root"
// common options
// --------------
//! \def SPC_CF_A_MAX_ITER
//! \brief Maximum iterations
#define SPC_CF_A_MAX_ITER "maximum_iterations"
//! \def SPC_CF_A_RES_TOL
//! \brief The residuals tolerance
//!
//! Absolute in SpecMiCP
//! Relative for DFPM
#define SPC_CF_A_RES_TOL "residual_tolerance"
//! \def SPC_CF_A_ABS_TOL
//! \brief Absolute residuals tolerance
#define SPC_CF_A_ABS_TOL "absolute_tolerance"
//! \def SPC_CF_A_STEP_TOL
//! \brief Absolute tolerance for the update
#define SPC_CF_A_STEP_TOL "step_tolerance"
//! \def SPC_CF_A_MAX_STEP_LENGTH
//! \brief Maximum value for the update
#define SPC_CF_A_MAX_STEP_LENGTH "maximum_step_length"
//! \def SPC_CF_A_MAX_STEP_MAX_ITER
//! \brief Maximum number at iterations when update value is maximum
#define SPC_CF_A_MAX_STEP_MAX_ITER "maximum_step_maximum_iterations"
//! \def SPC_CF_A_FILEPATH
//! \brief Path to a file
#define SPC_CF_A_FILEPATH "file"
//! \def SPC_CF_A_COMPONENT
//! \brief A component
#define SPC_CF_A_COMPONENT "component"
//! \def SPC_CF_A_NODES
//! \brief one or several nodes
//!
//! example: "1,2,5-8" for nodes 1,2,5,6,7 and 8
#define SPC_CF_A_NODES "nodes"
//! \def SPC_CF_A_TYPE
//! \brief Used when a thing can have several values
#define SPC_CF_A_TYPE "type"
//! \def SPC_CF_V_DEFAULT
//! \brief To indicate that a thing should be the default type
#define SPC_CF_V_DEFAULT "default"
//! \def SPC_CF_V_PLUGIN
//! \brief To indicate that a thing should use a user provided plugin
#define SPC_CF_V_PLUGIN "plugin"
//! \def SPC_CF_A_PLUGIN_FILE
//! \brief Attribute used to provide the path to a plugin file
#define SPC_CF_A_PLUGIN_FILE "plugin_file" // filepath to a plugin
//! \def SPC_CF_A_PLUGIN_NAME
//! \brief Attribute used to provide the name of a plugin
#define SPC_CF_A_PLUGIN_NAME "plugin_name" // name of the plugin to use
// user variables
// --------------
//! \def SPC_CF_S_USRVARS
//! \brief Section containing user defined variables
//!
//! These variables are used to compute algebraic expressions at runtime
#define SPC_CF_S_USRVARS "vars"
//! \def SPC_CF_S_USRVARS_SS_INT
//! \brief Subsection containing user defined integral variables
//!
//! Format is name: value
//!
//! Value can be an algebraic expression of previously defined integer variables
#define SPC_CF_S_USRVARS_SS_INT "int"
//! \def SPC_CF_S_USRVARS_SS_FLOAT
//! \brief Subsection containing user defined floating-point variables
//!
//! Format is name: value
//!
//! Value can be an algebraic expression of previously defined variables
//! (integer and floating points)
#define SPC_CF_S_USRVARS_SS_FLOAT "float"
// plugin manager
// --------------
//! \def SPC_CF_S_PLUGINS
//! \brief Section describing plugin options
#define SPC_CF_S_PLUGINS "plugins"
//! \def SPC_CF_S_PLUGINS_A_SEARCHDIRS
//! \brief Directories where plugins can be found
#define SPC_CF_S_PLUGINS_A_SEARCHDIRS "dirs"
//! \def SPC_CF_S_PLUGINS_A_TOLOAD
//! \brief Plugins to load
#define SPC_CF_S_PLUGINS_A_TOLOAD "to_load"
// Logs
// ----
// type of output
//! \def SPC_CF_S_LOGS_V_COUT
//! \brief Logs should be written to standard output
#define SPC_CF_S_LOGS_V_COUT "cout"
//! \def SPC_CF_S_LOGS_V_CERR
//! \brief Logs should be written to standard error output
#define SPC_CF_S_LOGS_V_CERR "cerr"
//! \def SPC_CF_S_LOGS_V_FILE
//! \brief Logs should be written to a file
#define SPC_CF_S_LOGS_V_FILE "file"
// the different levels
//! \def SPC_CF_S_LOGS_V_CRITICAL
//! \brief Log message should be at least of level "critical" to be written
#define SPC_CF_S_LOGS_V_CRITICAL "critical"
//! \def SPC_CF_S_LOGS_V_ERROR
//! \brief Log message should be at least of level "error" to be written
#define SPC_CF_S_LOGS_V_ERROR "error"
//! \def SPC_CF_S_LOGS_V_WARNING
//! \brief Log message should be at least of level "warning" to be written
#define SPC_CF_S_LOGS_V_WARNING "warning"
//! \def SPC_CF_S_LOGS_V_INFO
//! \brief Log message should be at least of level "info" to be written
#define SPC_CF_S_LOGS_V_INFO "info"
//! \def SPC_CF_S_LOGS_V_DEBUG
//! \brief Log message should be at least of level "debug" to be written
#define SPC_CF_S_LOGS_V_DEBUG "debug"
// runtime logs
//! \def SPC_CF_S_LOGS
//! \brief Runtime logs
#define SPC_CF_S_LOGS "logs"
//! \def SPC_CF_S_LOGS_A_LEVEL
//! \brief Output level for the runtime logs
#define SPC_CF_S_LOGS_A_LEVEL "level"
//! \def SPC_CF_S_LOGS_A_OUTPUT
//! \brief The ouput of the runtime logs
//!
//! Values are given by SPC_CF_S_LOGS_V_X where x is
//! either COUT, CERR, or FILE
#define SPC_CF_S_LOGS_A_OUTPUT "output"
//! \def SPC_CF_S_LOGS_A_FILEPATH
//! \brief the path to a file where logs should be written
//!
//! Only needed if logs output is a file
#define SPC_CF_S_LOGS_A_FILEPATH SPC_CF_A_FILEPATH
// configuration logs
//! \def SPC_CF_S_CONF_LOGS
//! \brief Logs written when configuration is parsed
#define SPC_CF_S_CONF_LOGS "configuration_logs"
//! \def SPC_CF_S_CONF_LOGS_A_OUTPUT
//! \brief The ouput of the configuration logs
//!
//! Values are given by SPC_CF_S_LOGS_V_X where x is
//! either COUT, CERR, or FILE
#define SPC_CF_S_CONF_LOGS_A_OUTPUT "output"
//! \def SPC_CF_S_CONF_LOGS_A_FILE
//! \brief the path to a file where logs configuration should be written
//!
//! Only needed if logs output is a file
#define SPC_CF_S_CONF_LOGS_A_FILE SPC_CF_A_FILEPATH
// units
// -----
//! \def SPC_CF_S_UNITS
//! \brief Configuration section for the units
#define SPC_CF_S_UNITS "units"
//! \def SPC_CF_S_UNITS_A_LENGTH
//! \brief Configure the length unit
#define SPC_CF_S_UNITS_A_LENGTH "length"
//! \def SPC_CF_S_UNITS_A_MASS
//! \brief Configure the mass unit
#define SPC_CF_S_UNITS_A_MASS "mass"
//! \def SPC_CF_S_UNITS_A_QUANTITY
//! \brief Configure the amount of substance unit
#define SPC_CF_S_UNITS_A_QUANTITY "quantity"
// mesh
// ----
//! \def SPC_CF_S_MESH
//! \brief Section to configure the mesh
#define SPC_CF_S_MESH "mesh"
//! \def SPC_CF_S_MESH_A_TYPE
//! \brief Type of the mesh
#define SPC_CF_S_MESH_A_TYPE "type"
//! \def SPC_CF_S_MESH_SS_UNIFMESH
//! \brief Subsection to configure a uniform mesh
#define SPC_CF_S_MESH_SS_UNIFMESH "uniform_mesh"
//! \def SPC_CF_S_MESH_SS_UNIFMESH_A_DX
//! \brief The length of an element
#define SPC_CF_S_MESH_SS_UNIFMESH_A_DX "dx"
//! \def SPC_CF_S_MESH_SS_UNIFMESH_A_NBNODES
//! \brief The number of nodes
#define SPC_CF_S_MESH_SS_UNIFMESH_A_NBNODES "nb_nodes"
//! \def SPC_CF_S_MESH_SS_UNIFMESH_A_SECTION
//! \brief Cross section of the sample
#define SPC_CF_S_MESH_SS_UNIFMESH_A_SECTION "section"
//! \def SPC_CF_S_MESH_SS_RAMPMESH
//! \brief Subsection to configure a ramp mesh
//!
//! Ramp mesh start by an increasing element size and finish with a plateau
//! at a max element size
#define SPC_CF_S_MESH_SS_RAMPMESH "ramp_mesh"
//! \def SPC_CF_S_MESH_SS_RAMPMESH_A_MIN_DX
//! \brief The minimum size of an element
#define SPC_CF_S_MESH_SS_RAMPMESH_A_MIN_DX "min_dx"
//! \def SPC_CF_S_MESH_SS_RAMPMESH_A_MAX_DX
//! \brief The maximum size of an element
#define SPC_CF_S_MESH_SS_RAMPMESH_A_MAX_DX "max_dx"
//! \def SPC_CF_S_MESH_SS_RAMPMESH_A_RAMP_LENGTH
//! \brief Total length of the ramp
#define SPC_CF_S_MESH_SS_RAMPMESH_A_RAMP_LENGTH "length_ramp"
//! \def SPC_CF_S_MESH_SS_RAMPMESH_A_PLATEAU_LENGTH
//! \brief Total lenght of the plateau
#define SPC_CF_S_MESH_SS_RAMPMESH_A_PLATEAU_LENGTH "length_plateau"
//! \def SPC_CF_S_MESH_SS_RAMPMESH_A_SECTION
//! \brief Cross section of the sample
#define SPC_CF_S_MESH_SS_RAMPMESH_A_SECTION "section"
//! \def SPC_CF_S_MESH_SS_UNIFAXISYMESH
//! \brief An axysimetric mesh, uniform (in radius, no volume)
//!
//! Either the number of nodes or the step can be provided.
#define SPC_CF_S_MESH_SS_UNIFAXISYMESH "uniform_axisymmetric"
//! \def SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_RADIUS
//! \brief Total radius of the sample
#define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_RADIUS "radius"
//! \def SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_NBNODES
//! \brief Number of nodes
#define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_NBNODES "nb_nodes"
//! \def SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_HEIGHT
//! \brief Height of the sample
#define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_HEIGHT "height"
//! \def SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_DX
//! \brief Lenght between two nodes
#define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_DX "dx"
// database
// --------
//! \def SPC_CF_S_DATABASE
//! \brief Section to configure the database
#define SPC_CF_S_DATABASE "database"
//! \def SPC_CF_S_DATABASE_A_PATH
//! \brief Path to teh database file
//!
//! Can be name of a file if search directories has been correctly set.
//! In particular, the database should search in directories provided by
//! the environment variable SPECMICP_DATABASE_PATH
#define SPC_CF_S_DATABASE_A_PATH "path"
//! \def SPC_CF_S_DATABASE_A_CHECKCOMPO
//! \brief If true, database will run a consistency check
//! (formula, electroneutrality...)
//!
//! It is overwritten by a flag in database file if provided.
#define SPC_CF_S_DATABASE_A_CHECKCOMPO "check_compo"
//! \def SPC_CF_S_DATABASE_A_SWAP
//! \brief Map containing species to be swapped
#define SPC_CF_S_DATABASE_A_SWAP "swap_components"
//! \def SPC_CF_S_DATABASE_A_SWAP_K_IN
//! \brief Aqueous species to add as a component
#define SPC_CF_S_DATABASE_A_SWAP_K_IN "in"
//! \def SPC_CF_S_DATABASE_A_SWAP_K_OUT
//! \brief Current compoment to remove from basis
#define SPC_CF_S_DATABASE_A_SWAP_K_OUT "out"
//! \def SPC_CF_S_DATABASE_A_REMOVE_GAS
//! \brief A list of gas to remove from the database
#define SPC_CF_S_DATABASE_A_REMOVE_GAS "remove_gas"
//! \def SPC_CF_S_DATABASE_A_REMOVE_SOLIDS
//! \brief A list of solids to remove from the database
#define SPC_CF_S_DATABASE_A_REMOVE_SOLIDS "remove_solids"
//! \def SPC_CF_S_DATABASE_A_REMOVE_SORBED
//! \brief A list of sorbed species to remove from the database
#define SPC_CF_S_DATABASE_A_REMOVE_SSITES "remove_ssites"
//! \def SPC_CF_S_DATABASE_A_REMOVE_SORBED
//! \brief A list of sorbed species to remove from the database
#define SPC_CF_S_DATABASE_A_REMOVE_SORBED "remove_sorbed"
//! \def SPC_CF_S_DATABASE_A_EXTRA_GAS
//! \brief a YAML-format string to add gas to the database
#define SPC_CF_S_DATABASE_A_EXTRA_GAS "extra_gas"
//! \def SPC_CF_S_DATABASE_A_EXTRA_SOLIDS
//! \brief a YAML-format string to add solid phases to the database
#define SPC_CF_S_DATABASE_A_EXTRA_SOLIDS "extra_solids"
//! \def SPC_CF_S_DATABASE_A_EXTRA_SSITES
//! \brief a YAML-format string to add surface sites to the database
#define SPC_CF_S_DATABASE_A_EXTRA_SSITES "extra_ssites"
//! \def SPC_CF_S_DATABASE_A_EXTRA_SORBED
//! \brief a YAML-format string to add sorbed species to the database
#define SPC_CF_S_DATABASE_A_EXTRA_SORBED "extra_sorbed"
//! \def SPC_CF_S_DATABASE_A_REMOVE_HALF_CELLS
//! \brief If true remove all half cells reaction from the database
//!
//! half-cells reactions are reaction involving the electron component
#define SPC_CF_S_DATABASE_A_REMOVE_HALF_CELLS "remove_half_cells"
//! \def SPC_CF_S_DATABASE_A_LIST_SOLIDS_TOKEEP
//! \brief a list of solid phases to keep in the database
//!
//! All solid phases not listed will be removed from the database
#define SPC_CF_S_DATABASE_A_LIST_SOLIDS_TOKEEP "list_solids_to_keep"
// SpecMiCP options
// ----------------
//! \def SPC_CF_S_SPECMICP
//! \brief Section to configure specmicp solver options
//!
//! May have a different name if several option sets are needed (ex. ReactMiCP),
//! but attributes are the same.
#define SPC_CF_S_SPECMICP "specmicp_options"
//! \def SPC_CF_S_SPECMICP_A_MAX_ITER
//! \brief The maximum number of iterations allowed
#define SPC_CF_S_SPECMICP_A_MAX_ITER SPC_CF_A_MAX_ITER
//! \def SPC_CF_S_SPECMICP_A_RES_TOL
//! \brief The absolute residual tolerance
#define SPC_CF_S_SPECMICP_A_RES_TOL SPC_CF_A_RES_TOL
//! \def SPC_CF_S_SPECMICP_A_STEP_TOL
//! \brief The update tolerance
#define SPC_CF_S_SPECMICP_A_STEP_TOL SPC_CF_A_STEP_TOL
//! \def SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH
//! \brief The maximum update allowed
//!
//! This option is used to avoid divergence
#define SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH
//! \def SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER
//! \brief The maximum number of iterations at maximum update allowed
//!
//! This option is used to detect tolerances
#define SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER
//! \def SPC_CF_S_SPECMICP_A_SCALING
//! \brief If true solver automatically scale the linear problem
#define SPC_CF_S_SPECMICP_A_SCALING "enable_scaling"
//! \def SPC_CF_S_SPECMICP_A_NONMONOTONE
//! \brief If true, the solver uses a nonmonotone linesearch algorithm
#define SPC_CF_S_SPECMICP_A_NONMONOTONE "enable_nonmonotone_linesearch"
//! \def SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION
//! \brief Multiplication applied to the descent direction when used as update
//!
//! Must be positive.
#define SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION "factor_descent_direction"
//! \def SPC_CF_S_SPECMICP_A_COND_CHECK
//! \brief Upper threshold for the condition number of the matrix
//!
//! if =-1, the check is not run
#define SPC_CF_S_SPECMICP_A_COND_CHECK "threshold_condition_check"
//! \def SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH
//! \brief Threshold used to check if linesearch is stuck in a loop
#define SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH "threshold_cycling_linesearch"
//! \def SPC_CF_S_SPECMICP_A_NONIDEAL_TOL
//! \brief Relative tolerance for the non-ideality model solution
#define SPC_CF_S_SPECMICP_A_NONIDEAL_TOL "non_ideality_tolerance"
//! \def SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER
//! \brief The maximum number of iteration to solve the non-ideality model
#define SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER "non_ideality_" SPC_CF_A_MAX_ITER
//! \def SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC
//! \brief Cutoff to include an equation into the system
#define SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC "cutoff_total_concentration"
//! \def SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION
//! \brief If the solver fails, the component concentrations will be reset to this value
//!
//! Value must be the log10 of the concentration (e.g. -3 instead of 10^-3)
#define SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION "restart_concentration"
//! \def SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC
//! \brief If the solver fails, the volume fraction of water will be reset to this value
#define SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC "restart_water_volume_fraction"
//! \def SPC_CF_S_SPECMICP_A_UNDER_RELAXATION
//! \brief Under relaxation factor for water volume fraction update
#define SPC_CF_S_SPECMICP_A_UNDER_RELAXATION "under_relaxation"
// SpecMiCP formulation
// --------------------
//! \def SPC_CF_S_FORMULATION
//! \brief Create a new chemical system
#define SPC_CF_S_FORMULATION "formulation"
//! \def SPC_CF_S_FORMULATION_A_SOLUTION
//! \brief Configure the water in the system
//!
//! A map expecting keys amount and unit
#define SPC_CF_S_FORMULATION_A_SOLUTION "solution"
//! \def SPC_CF_S_FORMULATION_A_AQUEOUS
//! \brief Configure the aqueous solution in the system
//!
//! A list of maps expecting keys amount, unit and label
#define SPC_CF_S_FORMULATION_A_AQUEOUS "aqueous"
//! \def SPC_CF_S_FORMULATION_A_MINERALS
//! \brief Configure the solid phases in the system
//!
//! A list of maps expecting keys amount, unit and label
#define SPC_CF_S_FORMULATION_A_MINERALS "solid_phases"
//! \def SPC_CF_S_FORMULATION_A_AMOUNT
//! \brief Provide the amount of a species to be added
//!
//! Meaning depends on the unit
#define SPC_CF_S_FORMULATION_A_AMOUNT "amount"
//! \def SPC_CF_S_FORMULATION_A_UNIT
//! \brief Provide the unit
//!
//! Available units depend on the type of the species being added
#define SPC_CF_S_FORMULATION_A_UNIT "unit"
//! \def SPC_CF_S_FORMULATION_A_LABEL
//! \brief Specify the species to be added
//!
//! A list of maps expecting keys amount, unit and label
#define SPC_CF_S_FORMULATION_A_LABEL "label"
//! \def SPC_CF_S_INITIALIZATION
//! \brief Initialization of AdimSmartSolver
#define SPC_CF_S_INITIALIZATION "initialization"
//! \def SPC_CF_S_INITIALIZATION_A_SOLUTION
//! \brief The initial volume fraction
#define SPC_CF_S_INITIALIZATION_A_SOLUTION "solution"
//! \def SPC_CF_S_INITIALIZATION_A_MINERALS
//! \brief Initialization of solid phases volume fractions
//!
//! A map {solid_phase, volume_fraction} is expected
#define SPC_CF_S_INITIALIZATION_A_MINERALS "solid_phases"
//! \def SPC_CF_S_INITIALIZATION_A_AQUEOUS
//! \brief Initialization of aqueous species
//!
//! A map {aq_species, volume_fraction} is expected.
//! If aq_species is "all", values are applied to all components.
#define SPC_CF_S_INITIALIZATION_A_AQUEOUS "aqueous"
-
-
// SpecMiCP constraints
// --------------------
//! \def SPC_CF_S_CONSTRAINTS
//! \brief Contraintes for the SpecMiCP solver
#define SPC_CF_S_CONSTRAINTS "constraints"
//! \def SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER
//! \brief A component that will be adjusted to keep the electroneutrality
#define SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER "charge_keeper"
//! \def SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY
//! \brief List of components with fixed activity
//!
//! List of maps {label, amount} or {label, amount_log}
#define SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY "fixed_activity"
//! \def SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY
//! \brief List of components with fixed molality
//!
//! List of maps {label, amount} or {label, amount_log}
#define SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY "fixed_molality"
//! \def SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY
//! \brief List of gas/components with fixed fugacity
//!
//! List of maps {label_gas, label_component, amount} or
//! {label_gas, label_component, amount_log}
#define SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY "fixed_fugacity"
//! \def SPC_CF_S_CONSTRAINTS_A_FIXEDSATURATION
//! \brief Set the system to have a fix saturation
//!
//! Expect the saturation (0<S<1)
#define SPC_CF_S_CONSTRAINTS_A_FIXEDSATURATION "fixed_saturation"
//! \def SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM
//! \brief If true, the system is saturated
#define SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM "saturated_system"
//! \def SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER
//! \brief If true, mass conservation for water is solved, the default
#define SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "conservation_water"
//! \def SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL
//! \brief The surface model
//!
//! Only equilibrium is implemented
#define SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL "surface_model"
//! \def SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION
//! \brief The concentration of surface sites
#define SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION "surface_site_concentration"
//! \def SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION
//! \brief The volume fraction of inert materials
#define SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION "inert_volume_fraction"
+//! \def SPC_CF_S_CONSTRAINTS_A_DISABLEIMMOBILE
+//! \brief Disable the immobile species in solver
+#define SPC_CF_S_CONSTRAINTS_A_DISABLEIMMOBILE "disable_immobile_species"
//! \def SPC_CF_S_CONSTRAINTS_A_LABEL
//! \brief Key in maps, label of a species
#define SPC_CF_S_CONSTRAINTS_A_LABEL "label"
//! \def SPC_CF_S_CONSTRAINTS_A_AMOUNT
//! \brief Key in maps, The amount of a species (in unit given by variables)
#define SPC_CF_S_CONSTRAINTS_A_AMOUNT "amount"
//! \def SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG
//! \brief Key in maps, The log10 amount of a species (in unit given by variables
#define SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "amount_log"
//! \def SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT
//! \brief Key in maps, Label of a component
#define SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT "label_component"
//! \def SPC_CF_S_CONSTRAINTS_A_LABEL_GAS
//! \brief Key in maps, Label of a gas
#define SPC_CF_S_CONSTRAINTS_A_LABEL_GAS "label_gas"
// Variables
// ---------
//! \def SPC_CF_S_REACTMICPVARIABLES
//! \brief Initialize reactmicp variables
//!
//! Usually done through a plugin
#define SPC_CF_S_REACTMICPVARIABLES "variables"
//! \def SPC_CF_S_REACTMICPVARIABLES_A_TYPE
//! \brief The method to use to initialize variables.
//!
//! If a plugin is used, it also needs the filepath and the name attributes.
#define SPC_CF_S_REACTMICPVARIABLES_A_TYPE SPC_CF_A_TYPE
// Boundary conditions
// -------------------
//! \def SPC_CF_S_BOUNDARYCONDITIONS
//! \brief The boundary conditions
//!
//! Details are dependant on context
//! e.g. for reactmicp unsaturation look in reactmicp/io/configuration_unsaturated
#define SPC_CF_S_BOUNDARYCONDITIONS "boundary_conditions"
// transport options
// -----------------
//! \def SPC_CF_S_DFPM
//! \brief Configuration for the dfpm solver
#define SPC_CF_S_DFPM "transport_options"
//! \def SPC_CF_S_DFPM_A_MAX_ITER
//! \brief The maximum number of iterations allowed
#define SPC_CF_S_DFPM_A_MAX_ITER SPC_CF_A_MAX_ITER
//! \def SPC_CF_S_DFPM_A_RES_TOL
//! \brief the relative residuals tolerance
#define SPC_CF_S_DFPM_A_RES_TOL SPC_CF_A_RES_TOL
//! \def SPC_CF_S_DFPM_A_ABS_TOL
//! \brief the absolute residuals tolerance
#define SPC_CF_S_DFPM_A_ABS_TOL SPC_CF_A_ABS_TOL
//! \def SPC_CF_S_DFPM_A_STEP_TOL
//! \brief The tolerance for the update
#define SPC_CF_S_DFPM_A_STEP_TOL SPC_CF_A_STEP_TOL
//! \def SPC_CF_S_DFPM_A_TRSHOLD_STATIONARY
//! \brief Detection of a stationnary point
//!
//! This condition is checked if update is less than step_tol but the residuals
//! is greater that the residuals tolerance
#define SPC_CF_S_DFPM_A_TRSHOLD_STATIONARY "threshold_stationary"
//! \def SPC_CF_S_DFPM_A_MAX_STEP_LENGTH
//! \brief The maximum update allowed
#define SPC_CF_S_DFPM_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH
//! \def SPC_CF_S_DFPM_A_MAX_STEP_MAX_ITER
//! \brief The maximum number of iterations allowed to have the maximum update
//!
//! Used to detect divergence
#define SPC_CF_S_DFPM_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER
//! \def SPC_CF_S_DFPM_A_SPARSE_SOLVER
//! \brief The sparse solver to use
#define SPC_CF_S_DFPM_A_SPARSE_SOLVER "sparse_solver"
//! \def SPC_CF_S_DFPM_A_SPARSE_PIVOTS_TRSHOLD
//! \brief Pivots threshold, usage depends on sparse solver used
#define SPC_CF_S_DFPM_A_SPARSE_PIVOTS_TRSHOLD "sparse_solver_pivots_threshold"
//! \def SPC_CF_S_DFPM_A_LINESEARCH
//! \brief The linesearch to use
#define SPC_CF_S_DFPM_A_LINESEARCH "linesearch"
//! \def SPC_CF_S_DFPM_A_QUASI_NEWTON
//! \brief Expect an integer, the number of iterations without updating the jacobian
#define SPC_CF_S_DFPM_A_QUASI_NEWTON "quasi_newton"
// ReactMiCP coupling solver
// -------------------------
//! \def SPC_CF_S_REACTMICP
//! \brief Reactmicp coupling algorithm options
#define SPC_CF_S_REACTMICP "reactmicp_options"
//! \def SPC_CF_S_REACTMICP_A_RES_TOL
//! \brief The relative residual tolerance
#define SPC_CF_S_REACTMICP_A_RES_TOL SPC_CF_A_RES_TOL
//! \def SPC_CF_S_REACTMICP_A_ABS_TOL
//! \brief The absolute residual tolerance
#define SPC_CF_S_REACTMICP_A_ABS_TOL SPC_CF_A_ABS_TOL
//! \def SPC_CF_S_REACTMICP_A_STEP_TOL
//! \brief The update tolerance
#define SPC_CF_S_REACTMICP_A_STEP_TOL SPC_CF_A_STEP_TOL
//! \def SPC_CF_S_REACTMICP_A_GOOD_ENOUGH_TOL
//! \brief Relaxed tolerance when max. iterations are reached
#define SPC_CF_S_REACTMICP_A_GOOD_ENOUGH_TOL "good_enough_tolerance"
//! \def SPC_CF_S_REACTMICP_A_MAX_ITER
//! \brief Maximum number of iterations allowed
#define SPC_CF_S_REACTMICP_A_MAX_ITER SPC_CF_A_MAX_ITER
//! \def SPC_CF_S_REACTMICP_A_IMPL_UPSCALING
//! \brief if true, includes upscaling stagger in SIA
#define SPC_CF_S_REACTMICP_A_IMPL_UPSCALING "implicit_upscaling"
// ReactMiCP output
// ----------------
//! \def SPC_CF_S_REACTOUTPUT
//! \brief Configuration of the main output file for Reactmicp
#define SPC_CF_S_REACTOUTPUT "reactmicp_output"
//! \def SPC_CF_S_REACTOUTPUT_A_FILEPATH
//! \brief Path to the output file
#define SPC_CF_S_REACTOUTPUT_A_FILEPATH SPC_CF_A_FILEPATH
//! \def SPC_CF_S_REACTOUTPUT_A_TYPE
//! \brief type of the output file
#define SPC_CF_S_REACTOUTPUT_A_TYPE "type"
//! \def SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5
//! \brief Output file is HDF5
#define SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5 "hdf5"
//! \def SPC_CF_S_REACTOUTPUT_A_TYPE_V_CSV
//! \brief output file is CSV
#define SPC_CF_S_REACTOUTPUT_A_TYPE_V_CSV "csv"
//! \def SPC_CF_S_REACTOUTPUT_A_POROSITY
//! \brief ReactMiCP saturated - output prosity
#define SPC_CF_S_REACTOUTPUT_A_POROSITY "porosity"
//! \def SPC_CF_S_REACTOUTPUT_A_PH
//! \brief ReactMiCP saturated - output pH
#define SPC_CF_S_REACTOUTPUT_A_PH "ph"
//! \def SPC_CF_S_REACTOUTPUT_A_DIFFUSIVITY
//! \brief ReactMiCP saturated - output diffusion coefficient
#define SPC_CF_S_REACTOUTPUT_A_DIFFUSIVITY "diffusivity"
//! \def SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL
//! \brief ReactMiCP saturated - ouput volume fraction on solid phases
#define SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL "volume_fraction_solid"
//! \def SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC
//! \brief ReactMiCP saturated - ouput total aqueous concentrations
#define SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC "total_aqueous_concentration"
//! \def SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC
//! \brief ReactMiCP saturated - ouput total solid concentrations
#define SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC "total_solid_concentration"
//! \def SPC_CF_S_REACTOUTPUT_A_TOT_CONC
//! \brief ReactMiCP saturated - ouput total concentrations
#define SPC_CF_S_REACTOUTPUT_A_TOT_CONC "total_concentration"
//! \def SPC_CF_S_REACTOUTPUT_A_SATINDEX
//! \brief ReactMiCP saturated - ouput volume fraction on solid phases
#define SPC_CF_S_REACTOUTPUT_A_SATINDEX "saturation_index"
//! \def SPC_CF_S_REACTOUTPUT_A_DATABASE
//! \brief Path to working database save file
#define SPC_CF_S_REACTOUTPUT_A_DATABASE "database"
// ReactMiCP simulation stuff
// --------------------------
//! \def SPC_CF_S_SIMULINFO
//! \brief Configure some information about simulation
#define SPC_CF_S_SIMULINFO "simulation"
//! \def SPC_CF_S_SIMULINFO_A_NAME
//! \brief Name of the simulation
#define SPC_CF_S_SIMULINFO_A_NAME "name"
//! \def SPC_CF_S_SIMULINFO_A_OUTPUTPREFIX
//! \brief Prefix to apply to solver output files
#define SPC_CF_S_SIMULINFO_A_OUTPUTPREFIX "output_prefix"
//! \def SPC_CF_S_SIMULINFO_A_PRINTITER
//! \brief If true, output debug information about every timesteps
//!
//! Can be quite verbose
#define SPC_CF_S_SIMULINFO_A_PRINTITER "print_iter_info"
//! \def SPC_CF_S_SIMULINFO_A_OUTPUTSTEP
//! \brief Output will be produced every 'output_step' seconds of simulation
#define SPC_CF_S_SIMULINFO_A_OUTPUTSTEP "output_step"
// ReactMiCP adaptive timestep
// ---------------------------
//! \def SPC_CF_S_TIMESTEPPER
//! \brief Configure the adaptative timestepper algorithm
#define SPC_CF_S_TIMESTEPPER "timestepper"
//! \def SPC_CF_S_TIMESTEP_A_LOWER_BOUND
//! \brief The minimum timestep allowed
#define SPC_CF_S_TIMESTEP_A_LOWER_BOUND "minimum_dt"
//! \def SPC_CF_S_TIMESTEP_A_UPPER_BOUND
//! \brief The maximum timestep allowed
#define SPC_CF_S_TIMESTEP_A_UPPER_BOUND "maximum_dt"
//! \def SPC_CF_S_TIMESTEP_A_RESTART_DT
//! \brief Timestep used after a critical failure
#define SPC_CF_S_TIMESTEP_A_RESTART_DT "restart_dt"
//! \def SPC_CF_S_TIMESTEP_A_LOWER_TARGET
//! \brief The lower bound target for the number of iterations
#define SPC_CF_S_TIMESTEP_A_LOWER_TARGET "lower_iterations_target"
//! \def SPC_CF_S_TIMESTEP_A_UPPER_TARGET
//! \brief The upper bound target for the number of iterations
#define SPC_CF_S_TIMESTEP_A_UPPER_TARGET "upper_iterations_target"
//! \def SPC_CF_S_TIMESTEP_A_AVERAGE_PARAMETER
//! \brief Parameter for the exponential moving average
#define SPC_CF_S_TIMESTEP_A_AVERAGE_PARAMETER "average_parameter"
//! \def SPC_CF_S_TIMESTEP_A_FACTOR_IF_FAILURE
//! \brief Multiplication factor for timestep after a failure
#define SPC_CF_S_TIMESTEP_A_FACTOR_IF_FAILURE "decrease_factor_if_failure"
//! \def SPC_CF_S_TIMESTEP_A_FACTOR_IF_MINIMUM
//! \brief Multiplication factor if error is minimized
#define SPC_CF_S_TIMESTEP_A_FACTOR_IF_MINIMUM "increase_factor_if_minumum"
//! \def SPC_CF_S_TIMESTEP_A_FACTOR_IF_DECREASE
//! \brief Multiplication factor if upper bound is reached
#define SPC_CF_S_TIMESTEP_A_FACTOR_IF_DECREASE "decrease_factor"
//! \def SPC_CF_S_TIMESTEP_A_FACTOR_IF_INCREASE
//! \brief Multiplication factor if lower bound is reached
#define SPC_CF_S_TIMESTEP_A_FACTOR_IF_INCREASE "increase_factor"
// ReactMiCP Staggers
// ------------------
//! \def SPC_CF_S_STAGGERS
//! \brief Configuration for ReactMiCP staggers
#define SPC_CF_S_STAGGERS "staggers"
//! \def SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER
//! \brief Configuration for ReactMiCP upscaling stagger
#define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER "upscaling"
//! \def SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE
//! \brief the type of upscaling stagger
//!
//! Usually a plugin
#define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE SPC_CF_A_TYPE
//! \def SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER
//! \brief Configuration section for the chemistry stagger
#define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER "chemistry"
//! \def SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE
//! \brief The type of chemistry stagger
#define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE SPC_CF_A_TYPE
//! \def SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM
//! \brief The chemistry stagger only solves equilibrium
#define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM "equilibrium"
//! \def SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS
//! \brief Options for the equilibrium chemistry stagger
#define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS "equilibrium_options"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER
//! \brief Configuration for ReactMiCP transport stagger
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER "transport"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION
//! \brief Merge saturation and water pressure equation (deprecated)
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION "merge_saturation"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS
//! \brief MERGE liquid and gas transport equations for aq. component (deprecated)
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS "merge_aqueous"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_CUTOFFRESIDUAL
//! \brief Cutoff for the first residual
//!
//! If less than the threshold, it is set to 1.
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_CUTOFFRESIDUAL "cutoff_R_0_squared"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS
//! \brief Default options for the dfpm solver
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS "default_options"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS
//! \brief Options for the saturation equations solver
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS "saturation_options"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS
//! \brief Options for the aqueous component liquid transport equations
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS "aqueous_options"
//! \def SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS
//! \brief Options for the gas diffusion equations
#define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS "gas_options"
// Run !
// -----
//! \def SPC_CF_S_RUN
//! \brief Configuration of simulation run
#define SPC_CF_S_RUN "run"
//! \def SPC_CF_S_RUN_A_RUNUNTIL
//! \brief run simulation until t='rununtil' s
#define SPC_CF_S_RUN_A_RUNUNTIL "run_until"
#endif // SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H
diff --git a/tests/specmicp/adim/adimensional_system_configuration.cpp b/tests/specmicp/adim/adimensional_system_configuration.cpp
index 2744a4b..9bd57a2 100644
--- a/tests/specmicp/adim/adimensional_system_configuration.cpp
+++ b/tests/specmicp/adim/adimensional_system_configuration.cpp
@@ -1,192 +1,260 @@
#include "catch.hpp"
#include "specmicp_database/database.hpp"
#include "specmicp_common/io/safe_config.hpp"
#include "specmicp_common/io/config_yaml_sections.h"
#include "specmicp_database/io/configuration.hpp"
#include "specmicp_common/physics/io/configuration.hpp"
#include "specmicp/io/configuration.hpp"
#include "specmicp/adimensional/adimensional_system_structs.hpp"
#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
TEST_CASE("Configuration", "[SpecMiCP],[io],[configuration]")
{
auto conf = specmicp::io::YAMLConfigFile::load("specmicp_conf_test.yaml");
specmicp::RawDatabasePtr raw_db = specmicp::io::configure_database(
conf.get_section(SPC_CF_S_DATABASE), {});
auto units_set = specmicp::io::configure_units(conf.get_section(SPC_CF_S_UNITS));
SECTION("Options")
{
specmicp::AdimensionalSystemSolverOptions opts;
specmicp::io::configure_specmicp_options(opts, units_set,
conf.get_section(SPC_CF_S_SPECMICP));
CHECK(opts.solver_options.fvectol == 1e-8);
CHECK(opts.solver_options.steptol == 1e-14);
CHECK(opts.solver_options.maxstep == 20.0);
CHECK(opts.system_options.non_ideality_tolerance == 1e-14);
}
SECTION("ReactantBox")
{
specmicp::ReactantBox reactant_box =
specmicp::io::configure_specmicp_reactant_box(
raw_db,
units_set,
conf.get_section("speciation").get_section(0)
);
auto constraints = reactant_box.get_constraints(true);
CHECK(raw_db->get_id_component_from_element("K") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("Na") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("Cl") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("Ca") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("Al") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("Si") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("C") != specmicp::no_species);
CHECK(raw_db->get_id_component_from_element("S") == specmicp::no_species);
CHECK(constraints.charge_keeper == raw_db->get_id_component("HO[-]"));
CHECK(constraints.inert_volume_fraction == 0.0);
CHECK(constraints.water_equation == specmicp::WaterEquationType::SaturatedSystem);
CHECK(constraints.fixed_molality_cs[0].id_component == raw_db->get_id_component("Cl[-]"));
CHECK(constraints.fixed_molality_cs[0].log_value == Approx(std::log10(0.2)));
}
SECTION("Solving from configuration")
{
specmicp::AdimensionalSystemSolverOptions opts;
specmicp::io::configure_specmicp_options(opts, units_set,
conf.get_section(SPC_CF_S_SPECMICP));
specmicp::ReactantBox reactant_box =
specmicp::io::configure_specmicp_reactant_box(
raw_db,
units_set,
conf.get_section("speciation").get_section(0)
);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_db, constraints, opts);
specmicp::Vector x;
solver.initialize_variables(x, 0.5, -6);
auto perf = solver.solve(x);
CHECK(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success);
}
}
auto conf_fixed_sat_str = R"plop(
database:
path: ../../data/cemdata.yaml
swap_components:
- out: H[+]
in: HO[-]
# - out: HCO3[-]
# in: CO2
- out: Al[3+]
in: Al(OH)4[-]
remove_gas: true
remove_half_cells: true
units:
length: decimeter
speciation:
- name: "Ca / CO2 system"
formulation:
solution:
amount: 0.8
unit: kg
aqueous:
- label: H2CO3
amount: 5.0
unit: mol/kg
- label: KOH
amount: 0.1
unit: mol/kg
- label: NaOH
amount: 0.1
unit: g
solid_phases:
- label: C3S
amount: 3.0
unit: mol/dm^3
- label: C2S
amount: 2.0
unit: mol/dm^3
- label: C3A
amount: 0.3
unit: mol/dm^3
constraints:
charge_keeper: "HO[-]"
fixed_molality:
- label: "Cl[-]"
amount: 0.2
fixed_saturation: 0.65
)plop";
TEST_CASE("Configuration - Fixed saturation", "[SpecMiCP],[io],[configuration]")
{
auto conf = specmicp::io::YAMLConfigFile::load_from_string(conf_fixed_sat_str);
specmicp::RawDatabasePtr raw_db = specmicp::io::configure_database(
conf.get_section(SPC_CF_S_DATABASE), {});
auto units_set = specmicp::io::configure_units(conf.get_section(SPC_CF_S_UNITS));
SECTION("ReactantBox")
{
specmicp::ReactantBox reactant_box =
specmicp::io::configure_specmicp_reactant_box(
raw_db,
units_set,
conf.get_section("speciation").get_section(0)
);
auto constraints = reactant_box.get_constraints(true);
CHECK(constraints.water_equation == specmicp::WaterEquationType::FixedSaturation);
CHECK(constraints.water_parameter == 0.65);
}
SECTION("Solving from configuration")
{
specmicp::AdimensionalSystemSolverOptions opts;
specmicp::ReactantBox reactant_box =
specmicp::io::configure_specmicp_reactant_box(
raw_db,
units_set,
conf.get_section("speciation").get_section(0)
);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_db, constraints, opts);
specmicp::Vector x;
solver.initialize_variables(x, 0.5, -6);
auto perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success);
auto sol = solver.get_raw_solution(x);
auto extr = specmicp::AdimensionalSystemSolutionExtractor(sol, raw_db, units_set);
REQUIRE(extr.saturation_water() == Approx(0.65).epsilon(1e-6));
}
}
+
+
+auto conf_immobile_species_str = R"plop(
+database:
+ path: ../../data/cemdata.yaml
+ swap_components:
+ - out: H[+]
+ in: HO[-]
+ - out: Al[3+]
+ in: Al(OH)4[-]
+ remove_gas: true
+ remove_half_cells: true
+
+units:
+ length: decimeter
+
+speciation:
+ - name: "test disable immobile species"
+ formulation:
+ solution:
+ amount: 0.5
+ unit: L
+ aqueous:
+ - label: Ca(OH)2
+ amount: 0.1
+ unit: mol/kg
+ - label: Na2SO4
+ amount: 0.5
+ unit: mol/kg
+ constraints:
+ charge_keeper: "HO[-]"
+ disable_immobile_species: true
+)plop";
+
+TEST_CASE("Configuration - disable immobile species", "[SpecMiCP],[io],[configuration]")
+{
+ auto conf = specmicp::io::YAMLConfigFile::load_from_string(conf_immobile_species_str);
+ specmicp::RawDatabasePtr raw_db = specmicp::io::configure_database(
+ conf.get_section(SPC_CF_S_DATABASE), {});
+
+ auto units_set = specmicp::io::configure_units(conf.get_section(SPC_CF_S_UNITS));
+
+ SECTION("Solving from configuration")
+ {
+ specmicp::AdimensionalSystemSolverOptions opts;
+ specmicp::ReactantBox reactant_box =
+ specmicp::io::configure_specmicp_reactant_box(
+ raw_db,
+ units_set,
+ conf.get_section("speciation").get_section(0)
+ );
+
+
+ auto constraints = reactant_box.get_constraints(true);
+ specmicp::AdimensionalSystemSolver solver(raw_db, constraints, opts);
+
+ specmicp::Vector x;
+ solver.initialize_variables(x, 0.5, -6);
+
+ auto perf = solver.solve(x);
+ REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success);
+ auto sol = solver.get_raw_solution(x);
+ auto extr = specmicp::AdimensionalSystemSolutionExtractor(sol, raw_db, units_set);
+ for (auto m: raw_db->range_mineral()) {
+ CHECK(extr.volume_fraction_mineral(m) == 0.0);
+ }
+ }
+}
diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp
index cdb005d..f4a107d 100644
--- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp
+++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp
@@ -1,675 +1,721 @@
#include "catch.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "specmicp/problem_solver/step_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_database/database.hpp"
#include "specmicp_common/physics/laws.hpp"
#include <iostream>
TEST_CASE("Reactant Box", "[units],[init],[formulation]")
{
SECTION("Units setting - SI")
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
auto raw_data= thedatabase.get_database();
auto units_set = specmicp::units::SI_units;
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_aqueous_species("CaCl2", 5.0, "kg");
reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg");
auto vector = reactant_box.get_total_concentration(false);
auto mass_sol = 0.5 * specmicp::laws::density_water(units_set);
CHECK(vector(0) == Approx( mass_sol
/ raw_data->molar_mass_basis(0, units_set)
));
auto id_ca = raw_data->get_id_component("Ca[2+]");
auto id_na = raw_data->get_id_component("Na[+]");
auto id_cl = raw_data->get_id_component("Cl[-]");
auto id_cacl2 = raw_data->get_id_compound("CaCl2");
CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2)));
CHECK(vector(id_na) == Approx(0.5*mass_sol));
CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca)));
}
SECTION("Units setting - mol L")
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
auto raw_data= thedatabase.get_database();
auto units_set = specmicp::units::SI_units;
units_set.length = specmicp::units::LengthUnit::decimeter;
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(0.5, "L");
reactant_box.add_aqueous_species("CaCl2", 5.0, "kg");
reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg");
auto vector = reactant_box.get_total_concentration(false);
auto mass_sol = 0.5 * specmicp::laws::density_water(units_set);
CHECK(vector(0) == Approx( mass_sol
/ raw_data->molar_mass_basis(0, units_set)
));
auto id_ca = raw_data->get_id_component("Ca[2+]");
auto id_na = raw_data->get_id_component("Na[+]");
auto id_cl = raw_data->get_id_component("Cl[-]");
auto id_cacl2 = raw_data->get_id_compound("CaCl2");
CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2)));
CHECK(vector(id_na) == Approx(0.5*mass_sol));
CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca)));
}
}
TEST_CASE("Solving adimensional system with problem setting", "[specmicp, MiCP, program, adimensional, solver]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
SECTION("Solving problem with dissolver - simpler problem") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
formulation.mass_solution = 0.8;
formulation.amount_minerals = {{"Portlandite", 10.0}};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
total_concentrations *= 1e3;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = raw_data->get_id_component("HO[-]");
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().solver_options.maxstep = 50.0;
solver.get_options().solver_options.maxiter_maxstep = 100;
solver.get_options().solver_options.use_crashing = false;
solver.get_options().solver_options.use_scaling = true;
solver.get_options().solver_options.factor_descent_condition = 1e-10;
solver.get_options().solver_options.factor_gradient_search_direction = 200;
specmicp::Vector x;
solver.initialize_variables(x, 0.8, -2);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set);
CHECK(extr.volume_fraction_water() == Approx(0.802367).epsilon(1e-4));
CHECK(extr.volume_fraction_mineral(raw_data->get_id_mineral("Portlandite")) == Approx(0.32947).epsilon(1e-2));
}
SECTION("Solving problem with dissolver") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
formulation.mass_solution = 0.8;
formulation.amount_minerals = {{"C3S", 0.7}, {"C2S", 0.3}};
formulation.extra_components_to_keep = {"HCO3[-]", };
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
total_concentrations *= 1e3;
total_concentrations(2) = 500;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = raw_data->get_id_component("HO[-]");
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().solver_options.maxstep = 20.0;
solver.get_options().solver_options.maxiter_maxstep = 100;
solver.get_options().solver_options.use_crashing = false;
solver.get_options().solver_options.use_scaling = true;
solver.get_options().solver_options.factor_descent_condition = 1e-10;
solver.get_options().solver_options.factor_gradient_search_direction = 200;
specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral());
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
}
SECTION("Solving problem with reactant box") {
std::cout << "\n Solving with reactant box \n ------------------------ \n";
std::cout << "== with SI units == \n";
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
reactant_box.set_charge_keeper("HO[-]");
specmicp::AdimensionalSystemConstraints constraints = reactant_box.get_constraints(true);
specmicp::Vector total_concentrations = constraints.total_concentrations;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.5, -6);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(solution,raw_data, specmicp::units::SI_units);
std::cout << "== with non SI units : mmol/dm == \n";
specmicp::units::UnitsSet dmmol;
dmmol.length = specmicp::units::LengthUnit::decimeter;
dmmol.quantity = specmicp::units::QuantityUnit::millimoles;
reactant_box.set_units(dmmol);
specmicp::AdimensionalSystemConstraints constraints2 = reactant_box.get_constraints(false);
specmicp::Vector total_concentrations2 = constraints2.total_concentrations;
specmicp::AdimensionalSystemSolver solver2(raw_data, constraints2);
solver2.get_options().units_set = dmmol;
solver2.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::Vector x2;
solver2.initialize_variables(x2, 0.5, -6);
perf = solver2.solve(x2);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
specmicp::AdimensionalSystemSolution solution2 = solver2.get_raw_solution(x);
CHECK(solution.main_variables(0) == Approx(solution2.main_variables(0)).epsilon(1e-6));
CHECK(solution.main_variables(2) == Approx(solution2.main_variables(2)).epsilon(1e-6));
CHECK(solution.main_variables(5) == Approx(solution2.main_variables(5)).epsilon(1e-6));
CHECK(solution.ionic_strength == Approx(solution2.ionic_strength).epsilon(1e-6));
specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, dmmol);
for (auto idc: raw_data->range_component()) {
if (idc == specmicp::database::electron_id) continue;
// std::cout << "Id " << idc << " - " << raw_data->get_label_component(idc) << std::endl;
// std::cout << total_concentrations(idc) << " - "
// << total_concentrations2(idc) << " - "
// << extr.total_concentration(idc) << " ~ "
// << extr.volume_fraction_gas_phase() * extr.total_gaseous_concentration(idc) << " # "
// << extr.fugacity_gas(0) << " - "
// << extr2.total_concentration(idc)
// << std::endl;
CHECK(extr.total_concentration(idc) == Approx(total_concentrations(idc)).epsilon(1e-8));
CHECK(extr2.total_concentration(idc) == Approx(total_concentrations2(idc)).epsilon(1e-8));
CHECK(extr.total_concentration(idc) == Approx(extr2.total_concentration(idc)).epsilon(1e-8));
}
std::cout << "== with non SI units : mmol/cm == \n";
specmicp::units::UnitsSet cmmol;
cmmol.length = specmicp::units::LengthUnit::centimeter;
cmmol.quantity = specmicp::units::QuantityUnit::millimoles;
reactant_box.set_units(cmmol);
specmicp::AdimensionalSystemConstraints constraints3 = reactant_box.get_constraints(false);
specmicp::Vector total_concentrations3 = constraints3.total_concentrations;
specmicp::AdimensionalSystemSolver solver3(raw_data, constraints3);
solver3.get_options().units_set = cmmol;
solver3.get_options().solver_options.set_tolerance(1e-8, 1e-12);
specmicp::Vector x3;
solver3.initialize_variables(x3, 0.5, -6);
perf = solver3.solve(x3);
for (auto idc: raw_data->range_component()) {
if (idc == specmicp::database::electron_id) continue;
CHECK(total_concentrations(idc) == Approx(1e3*total_concentrations3(idc)).epsilon(1e-8));
}
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
specmicp::AdimensionalSystemSolution solution3 = solver3.get_raw_solution(x);
CHECK(solution.main_variables(0) == Approx(solution3.main_variables(0)).epsilon(1e-6));
CHECK(solution.main_variables(2) == Approx(solution3.main_variables(2)).epsilon(1e-6));
CHECK(solution.main_variables(5) == Approx(solution3.main_variables(5)).epsilon(1e-6));
CHECK(solution.ionic_strength == Approx(solution3.ionic_strength).epsilon(1e-6));
specmicp::AdimensionalSystemSolutionExtractor extr3(solution3, raw_data, cmmol);
for (auto idc: raw_data->range_component()) {
CHECK(extr3.total_concentration(idc) == Approx(total_concentrations3(idc)).epsilon(1e-8));
}
}
+ SECTION("Solving problem with immobilized immobile species") {
+
+ specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
+ std::map<std::string, std::string> swapping ({
+ {"H[+]","HO[-]"},
+ {"Si(OH)4", "SiO(OH)3[-]"},
+ {"Al[3+]", "Al(OH)4[-]"},
+ {"HCO3[-]", "CO2"}
+ });
+ thedatabase.swap_components(swapping);
+ //thedatabase.remove_gas_phases();
+ thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
+ specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
+
+ specmicp::units::UnitsSet dm_mol;
+ dm_mol.length = specmicp::units::LengthUnit::decimeter;
+ specmicp::ReactantBox reactant_box(raw_data, dm_mol);
+
+ reactant_box.set_solution(0.5, "L");
+ reactant_box.add_aqueous_species("Ca(OH)2", 0.1, "mol/kg");
+ reactant_box.add_aqueous_species("Na2SO4", 0.5, "mol/kg");
+ reactant_box.disable_immobile_species();
+ reactant_box.set_charge_keeper("HO[-]");
+
+ specmicp::AdimensionalSystemConstraints constraints = reactant_box.get_constraints(true);
+ specmicp::Vector total_concentrations = constraints.total_concentrations;
+
+ specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
+
+ solver.get_options().units_set = dm_mol;
+ solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
+
+ specmicp::Vector x;
+ solver.initialize_variables(x, 0.5, -6);
+ auto perf = solver.solve(x);
+ REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
+
+ auto solution = solver.get_raw_solution(x);
+ specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, dm_mol);
+
+ for (auto m: raw_data->range_mineral()) {
+ CHECK(extr.volume_fraction_mineral(m) == 0.0);
+ CHECK(extr.saturation_index(m) > 0.0);
+ }
+ }
+
SECTION("Solving problem with reactant box and smart solver") {
std::cout << "\n Solving with reactant box and smart solver \n ------------------------ \n";
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
reactant_box.set_charge_keeper("HO[-]");
specmicp::SmartAdimSolver solver(raw_data, reactant_box);
solver.set_init_volfrac_water(0.5);
solver.set_init_molality(1e-6);
solver.solve();
REQUIRE(solver.is_solved() == true);
}
SECTION("Solving problem with reactant box and smart solver, fixed saturation") {
std::cout << "\n Solving with reactant box and smart solver, fixed saturation \n ------------------------ \n";
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
reactant_box.set_charge_keeper("HO[-]");
reactant_box.set_fixed_saturation(0.8);
specmicp::SmartAdimSolver solver(raw_data, reactant_box);
solver.set_init_volfrac_water(0.5);
solver.set_init_molality(1e-6);
solver.solve();
REQUIRE(solver.is_solved() == true);
auto sol = solver.get_solution();
auto extr = specmicp::AdimensionalSystemSolutionExtractor(sol, raw_data, reactant_box.get_units());
CHECK(extr.saturation_water() == Approx(0.8).epsilon(1e-6));
}
SECTION("Reactant box fixed activity") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(0.5, "L/L");
reactant_box.add_solid_phase("C3S", 0.350, "L/L");
reactant_box.add_solid_phase("C2S", 0.05, "L/L");
reactant_box.set_charge_keeper("HCO3[-]");
reactant_box.add_fixed_activity_component("H[+]", 1e-9);
auto constraints = reactant_box.get_constraints(true);
CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]"));
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.6, -4);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.pH() == Approx(9));
}
SECTION("Reactant box fixed molality") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"},
// {"H[+]", "HO[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(0.5, "L/L");
reactant_box.add_solid_phase("C3S", 0.35, "L/L");
reactant_box.add_solid_phase("C2S", 0.05, "L/L");
reactant_box.set_charge_keeper("HCO3[-]");
reactant_box.add_fixed_molality_component("H[+]", 1e-5);
auto constraints = reactant_box.get_constraints(true);
CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]"));
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.6, -4);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.molality_component(raw_data->get_id_component("H[+]")) == Approx(1e-5));
}
SECTION("Reactant box fixed fugacity") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_fixed_fugacity_gas("CO2(g)", "HCO3[-]", 1e-3);
reactant_box.set_charge_keeper("H[+]");
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.8, -6);
solver.get_options().solver_options.set_tolerance(1e-6, 1e-8);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.fugacity_gas(raw_data->get_id_gas("CO2(g)")) == Approx(1e-3));
}
SECTION("Reactant box - mineral upper bound") {
std::cout << "\n Reactant box : box vi \n ------------------------ \n";
std::cout << "first solution" << std::endl;
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
//{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]", "Al(OH)4[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
thedatabase.minerals_keep_only({
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2(am)",
"Calcite",
"Al(OH)3(mic)",
"C3AH6",
"C3ASO_41H5_18",
"C3ASO_84H4_32",
"C4AH19",
"C4AH13",
"C2AH7_5",
"CAH10",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Gypsum",
"Ettringite"
});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::units::UnitsSet units_set;
units_set.length = specmicp::units::LengthUnit::centimeter;
units_set.quantity = specmicp::units::QuantityUnit::millimoles;
specmicp::ReactantBox reactant_box(raw_data, units_set);
reactant_box.set_solution(0.25, "kg/L");
reactant_box.add_solid_phase("Portlandite", 4.46, "mol/L");
reactant_box.add_solid_phase("CSH,jennite", 4.59, "mol/L");
reactant_box.add_solid_phase("Monosulfoaluminate", 0.3, "mol/L");
reactant_box.add_solid_phase("Al(OH)3(mic)", 0.6, "mol/L");
reactant_box.add_solid_phase("Calcite", 0.2, "mol/L");
//reactant_box.add_aqueous_species("CO2", 0.5, "mol/L");
reactant_box.add_aqueous_species("KOH", 473.0, "mmol/L");
reactant_box.add_aqueous_species("NaOH", 52.0, "mmol/L");
reactant_box.set_fixed_saturation(0.7);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().units_set = units_set;
solver.get_options().system_options.non_ideality_tolerance = 1e-12;
specmicp::Vector x;
solver.initialize_variables(x, 0.7, -6);
// /solver.get_options().solver_options.set_tolerance(1e-6, 1e-8);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, units_set);
std::cout << "constrained solution" << std::endl;
specmicp::index_t id_ettr = raw_data->get_id_mineral("Ettringite");
specmicp::scalar_t max_vol_frac = extr.volume_fraction_mineral(id_ettr);
reactant_box.set_mineral_upper_bound("Ettringite", max_vol_frac);
reactant_box.add_solid_phase("Gypsum", 0.5, "mol/L");
constraints = reactant_box.get_constraints(false);
solver = specmicp::AdimensionalSystemSolver(raw_data, constraints);
solver.get_options().units_set = units_set;
perf = solver.solve(x, false);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol2 = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr2(sol2, raw_data, units_set);
CHECK(extr2.volume_fraction_mineral(id_ettr) <= max_vol_frac);
std::cout << extr2.volume_fraction_mineral(id_ettr) << " ?(<=) " << max_vol_frac << std::endl;
}
SECTION("Solving problem with stepper") {
std::cout << "\n Solving with stepper \n ------------------------ \n";
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
// initial constraints
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.set_charge_keeper("HO[-]");
reactant_box.add_component("CO2");
specmicp::uindex_t max_step = 20;
specmicp::StepProblem problem;
problem.set_initial_constraints(reactant_box, true);
specmicp::index_t id_co2 = raw_data->get_id_component("CO2");
specmicp::index_t id_ca = raw_data->get_id_component("Ca[2+]");
// The updates to apply at each step
specmicp::Vector update;
update.setZero(raw_data->nb_component());
update(id_co2) = (0.95*problem.get_constraints().total_concentrations(id_ca))/max_step;
problem.set_constraint_updater(specmicp::StepTotalConcentration(update));
problem.set_stop_condition_step(max_step);
specmicp::StepProblemRunner runner;
// output
std::vector<double> pHs;
pHs.reserve(max_step);
runner.set_step_output(
[&pHs](specmicp::uindex_t cnt,
const specmicp::AdimensionalSystemSolutionExtractor& extr) {
pHs.push_back(extr.pH());
return specmicp::StepProblemStatus::OK;
});
auto status = runner.run(problem);
REQUIRE(status == specmicp::StepProblemStatus::OK);
CHECK(pHs.size() == max_step);
for (auto it: pHs) {
std::cout << it << std::endl;
}
CHECK(pHs[pHs.size()-1] == Approx(9.84065).epsilon(1e-4));
}
}