diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index f7c7c2a..550a8cd 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,1933 +1,2043 @@ /* ============================================================================= Copyright (c) 2014-2017 F. Georget Princeton University Copyright (c) 2017-2018 F. Georget 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 #include -//#include +#include // 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(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), - m_second(ptrdata), + m_second(ptrdata.get()), 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(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), - m_second(previous_solution), + m_second(previous_solution, ptrdata.get()), 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 + database::DataContainer* 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())) + sorbed_concentrations(Vector::Zero(data->nb_sorbed())), + logactivity_solid_phases(Vector::Zero(data->nb_mineral())), + molar_fractions(Vector::Zero(data->nb_solid_solutions())) {} // Previous solution // ----------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( - const AdimensionalSystemSolution& previous_solution + const AdimensionalSystemSolution& previous_solution, + database::DataContainer* data ): 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) -{} + sorbed_concentrations(previous_solution.sorbed_molalities), + logactivity_solid_phases(Vector::Zero(data->nb_mineral())), + molar_fractions(Vector::Zero(data->nb_solid_solutions())) +{ + // ToDo : initialization of solid solutions +} // 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()) + active_sorbed(data->nb_sorbed()), + active_ssol(data->nb_solid_solutions(), true) {} // 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(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(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(SurfaceEquationType::NoEquation); } } // No surface potential equation for this model if (at_least_one) { m_equations.type_equation(dof_surface_potential()) = static_cast(SurfaceEquationType::Equilibrium); } else { m_equations.type_equation(dof_surface_potential()) = static_cast(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(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(SurfaceEquationType::NoEquation); } } if (at_least_one) { m_equations.type_equation(dof_surface_potential()) = static_cast(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(SurfaceEquationType::NoEquation); } } else { m_equations.type_equation(dof_surface_potential()) = static_cast(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 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(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; + // m_equations.nb_free_variables = neq; // following equations are complementarity conditions // Minerals // ======== - if (not constraints.immobile_disabled) - { + if (constraints.immobile_disabled) { + m_equations.nb_free_variables = neq; + } else { + std::vector can_precipitate(m_data->nb_mineral(), true); for (index_t m: m_data->range_mineral()) { - bool can_precipitate = true; // true if all the components of the mineral + //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) { if (m_data->nu_mineral(m, k) != 0.0 and molar_volume > 0.0) { - can_precipitate = false; + can_precipitate[m] = 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) { if (v.id_mineral == m) { if (v.param == 0) { - can_precipitate = false; + can_precipitate[m] = false; } } } } // Not an equation if fixed saturation index for (const auto& it : constraints.fixed_SI_cs) { if (it.id_mineral == m) { - can_precipitate = false; + can_precipitate[m] = false; + } + } + } + if (get_options().solve_solid_solutions) { + for (auto i: m_data->range_solid_solutions()) + { + for (auto m=m_data->members_solid_solution(i);m;++m) { + specmicp_assert(can_precipitate[m.index()] = true); + m_equations.add_equation(dof_mineral(m.index()), &neq); + ++m_equations.nb_free_variables; } } - if (can_precipitate) + } + // above equations are 'free' (i.e. non constrained) + m_equations.nb_free_variables = neq; + // following equations are complementarity conditions + for (index_t m: m_data->range_mineral()) + { + if (can_precipitate[m] and ideq_min(m)==no_equation) + // can precipitate and not part of solid solution { m_equations.add_equation(dof_mineral(m), &neq); } } + //if (can_precipitate) + //{ + // m_equations.add_equation(dof_mineral(m), &neq); + // } else if (m_data->nb_solid_solutions() > 0) { + // index_t ssol = m_data->id_solid_solution_of_mineral(m); + // m_equations.set_solid_solution_active(ssol, false); + //} + } 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.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(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(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(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(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(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(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(EqT::FixedMolality); m_fixed_values(it.id_component) = it.log_value; } for (const auto& it: constraints.fixed_SI_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '"+ m_data->components.get_label(it.id_component) + "' is already constraint, a fixed saturation index condition can not be applied."); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedSaturationIndex); m_fixed_values(it.id_component) = it.saturation_index; m_equations.related_species(it.id_component) = it.id_mineral; } // 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(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(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_fixed_saturation_index(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedSaturationIndex); const index_t id_m = m_equations.fixed_activity_species[component]; scalar_t res = fixed_saturation_index_bc(component) + m_data->logk_mineral(id_m); for (index_t component: m_data->range_aqueous_component()) { if (m_data->nu_mineral(id_m, component) == 0) continue; res -= m_data->nu_mineral(id_m, component)*( log_gamma_component(component) + log_component_molality(x, component) ); } // No scaling because Fixed Saturation Index is usually small and often 0.0 ! //res /= fixed_saturation_index_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); + scalar_t res = m_data->logk_mineral(m) + m_second.logactivity_solid_phases(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); + if (get_options().solve_solid_solutions) + compute_solid_solutions(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; case AqueousComponentEquationType::FixedSaturationIndex: residual(ideq_paq(i)) = residual_fixed_saturation_index(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; jmolar_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; } // Fixed saturation index // ====================== case AqueousComponentEquationType::FixedSaturationIndex: { const auto id_m = m_equations.fixed_activity_species[i]; for (auto k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation or m_data->nu_mineral(id_m, k) == 0.0) continue; jacobian(idp, ideq_paq(k)) = -m_data->nu_mineral(id_m, k); } 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()); + if (get_options().solve_solid_solutions) { + index_t ids = m_data->id_solid_solution_of_mineral(m); + if (ids != no_species) { + scalar_t tot_c = 0; + for (auto it=m_data->members_solid_solution(ids); it; ++it) { + tot_c += ( volume_fraction_mineral(x, it.index()) + /m_data->molar_volume_mineral(it.index()) ); + //std::cout << volume_fraction_mineral(x, it.index()) << " - "; + } + //std::cout << "p " << tot_c << std::endl; + if (tot_c > 0) { + for (auto it=m_data->members_solid_solution(ids); it; ++it) { + auto mc = it.index(); + auto idmc = ideq_min(mc); + //std::cout << m_data->get_label_mineral(mc) << " =-= " + // << mc << " -=- " << idmc << std::endl; + jacobian(idm, idmc) += 1.0/std::pow(10, m_second.logactivity_solid_phases(mc)) *( + tot_c- (volume_fraction_mineral(x, mc) + /m_data->molar_volume_mineral(mc)))/std::pow(tot_c, 2); + } + } + } + } + } } 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); + if (get_options().solve_solid_solutions) + compute_solid_solutions(x); +} + +void AdimensionalSystem::compute_solid_solutions(const Vector& x) +{ + m_second.logactivity_solid_phases.setZero(); + for (auto id: m_data->range_solid_solutions()) { + if (not m_equations.active_ssol[id]) continue; + scalar_t tot_c = 0; + for (auto it=m_data->members_solid_solution(id); it; ++it) { + tot_c += ( volume_fraction_mineral(x, it.index()) + /m_data->molar_volume_mineral(it.index()) ); + } + if (tot_c > 1e-10) { + for (auto it=m_data->members_solid_solution(id); it; ++it) { + auto m = it.index(); + m_second.logactivity_solid_phases(m) = std::log10( + (volume_fraction_mineral(x, m) + /m_data->molar_volume_mineral(m))/tot_c); + } + } + } } 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) { + 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) { + //std::cout << get_options().solve_solid_solutions << std::endl; 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 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); + if (get_options().solve_solid_solutions) + { + for (auto d: m_data->range_solid_solutions()) + { + for (auto m=m_data->members_solid_solution(d);m;++m) + { + xtot(dof_mineral(m)) = 0.1; + } + } + } + m_second = SecondaryVariables(m_data.get()); } 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); + if (get_options().solve_solid_solutions) + { + for (auto d: m_data->range_solid_solutions()) + { + for (auto m=m_data->members_solid_solution(d);m;++m) + { + xtot(dof_mineral(m)) = 0.1; + } + } + } + m_second = SecondaryVariables(m_data.get()); } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp index e790200..ba7d6bb 100644 --- a/src/specmicp/adimensional/adimensional_system.hpp +++ b/src/specmicp/adimensional/adimensional_system.hpp @@ -1,894 +1,907 @@ /* ============================================================================= Copyright (c) 2014-2017 F. Georget Princeton University Copyright (c) 2017-2018 F. Georget 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/maths.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, public AdimemsionalSystemNumbering, public OptionsHandler, 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 correspoinding to a fixed saturation index 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_saturation_index(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(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(m_equations.type_equation(dof_component(component))); } //! \brief Return the type of equation for the electron ElectronEquationType electron_equation_type() const { return static_cast(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(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& 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_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 fixed saturation index for 'component' scalar_t fixed_saturation_index_bc(index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedSaturationIndex); 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 Compute the solid solutions + //! + void compute_solid_solutions(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 + Vector molar_fractions; //!< The molar fractions of solid solutions + Vector logactivity_solid_phases; //!< The log of activity of solid phases + //! \brief Initialization without a previous solution - SecondaryVariables(const RawDatabasePtr& data); + SecondaryVariables(database::DataContainer* data); //! \brief Initialization with a previous solution - SecondaryVariables(const AdimensionalSystemSolution& previous_solution); + SecondaryVariables(const AdimensionalSystemSolution& previous_solution, + database::DataContainer* data); }; //! \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 ideq; std::vector component_equation_type; std::vector fixed_activity_species; std::vector nonactive_component {}; std::vector active_aqueous; std::vector active_gas; std::vector active_sorbed; + std::vector active_ssol; bool is_box_constrained {false}; std::vector is_box_vi {}; std::vector 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; } + //! \brief Set the active flag of solid solution 'id' to 'is_active' + void set_solid_solution_active(index_t id, bool is_active) { + active_ssol[id] = is_active; + } + //! \brief Initialise the box constrained data void activate_box_constrained(index_t nb_eq) { is_box_constrained = true; is_box_vi = std::vector(nb_eq, false); upper_bounds = std::vector(nb_eq, 0.0); } - + //! \brief Set the mineral equation to be a box constrained equation 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_numbering.hpp b/src/specmicp/adimensional/adimensional_system_numbering.hpp index 5c7b324..daa9081 100644 --- a/src/specmicp/adimensional/adimensional_system_numbering.hpp +++ b/src/specmicp/adimensional/adimensional_system_numbering.hpp @@ -1,134 +1,144 @@ /* ============================================================================= Copyright (c) 2014 - 2017 F. Georget Princeton University 2017 - 2018 F. Georget 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_ADIMENSIONALSYSTEMNUMBERING #define SPECMICP_ADIMENSIONALSYSTEMNUMBERING #include "specmicp_common/types.hpp" #ifndef SPECMICP_DATABASE_HPP #include "specmicp_database/database.hpp" #endif namespace specmicp { //! \brief The equation numbering scheme for the adimensional system //! //! This is the numbering of the main variables (the non-reduced vector) class AdimemsionalSystemNumbering { public: //! \brief Initialize the numbering scheme //! \param ptr_database shared_ptr to the database container AdimemsionalSystemNumbering(RawDatabasePtr ptr_database): m_data(ptr_database) {} //! \brief Return the dof for water conservation index_t dof_water() const noexcept {return m_data->water_index();} //! \brief Return the dof corresponding to the component conservation //! \param component Index of the component index_t dof_component(index_t component) const NOEXCEPT { specmicp_assert_component_bounds(component, m_data); return component; } //! \brief Return the dof corresponding to the electron activity index_t dof_electron() const noexcept { return m_data->electron_index(); } //! \brief Return the offset for the dof of sorbed species index_t offset_surface() const noexcept { return m_data->nb_component(); } //! \brief Return the dof corresponding to the sorption sites conservation index_t dof_surface(index_t q) const noexcept { return offset_surface()+q; } //! \brief Return the offset of the surface potential index_t offset_surface_potential() const noexcept { return offset_surface()+m_data->nb_ssites(); } //! \brief Return the doff for the surface potential index_t dof_surface_potential() const noexcept { return offset_surface_potential(); } //! \brief Return the offset for the dof of minerals //! \sa dof_mineral index_t offset_minerals() const noexcept { return offset_surface_potential()+1; } //! \brief Return the dof for mineral 'mineral' //! \param mineral Index of the mineral (in the database) index_t dof_mineral(index_t mineral) const NOEXCEPT { specmicp_assert_mineral_bounds(mineral, m_data); return offset_minerals()+mineral; } - + //! \brief Return the offset for the dof of solid solution + //! \sa dof_solid_solution + index_t offset_solid_solution() const noexcept { + return offset_minerals()+m_data->nb_mineral(); + } + //! \brief Return the dof for solid solution 'ssol' + index_t dof_solid_solution(index_t ssol) const NOEXCEPT { + specmicp_assert(ssol>0 and ssolnb_solid_solutions()); + return offset_solid_solution()+ssol; + } //! \brief Return the number of dofs index_t total_dofs() const noexcept { return offset_minerals()+m_data->nb_mineral(); + //return offset_solid_solution()+m_data->nb_solid_solutions(); } // secondary species //! \brief dof of the activitiy coefficient for 'component' //! \param component Index of the component (in the database) - index_t dof_component_gamma(index_t component) const { + index_t dof_component_gamma(index_t component) const NOEXCEPT { specmicp_assert_component_bounds(component, m_data); return component; } //! \brief dof of a secondary aqueous species //! \param aqueous Index of the aqueous species (in the database) - index_t dof_aqueous(index_t aqueous) const { + index_t dof_aqueous(index_t aqueous) const NOEXCEPT { specmicp_assert(aqueous < m_data->nb_aqueous()); return aqueous; } //! \brief dof of the activity coefficient for 'aqueous' //! \param aqueous Index of the aqueous species (in the database) - index_t dof_aqueous_gamma(index_t aqueous) const { + index_t dof_aqueous_gamma(index_t aqueous) const NOEXCEPT { specmicp_assert_aqueous_bounds(aqueous, m_data); return m_data->nb_component()+aqueous; } //! \brief Return the database RawDatabasePtr get_database() { return m_data; } protected: RawDatabasePtr m_data; //!< The database }; } // end namespace specmicp #endif // SPECMICP_ADIMENSIONALSYSTEMNUMBERING diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp index 3402e30..702afc4 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.cpp +++ b/src/specmicp/adimensional/adimensional_system_solver.cpp @@ -1,551 +1,594 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "adimensional_system_solver.hpp" #include "adimensional_system.hpp" #include "adimensional_system_solution.hpp" #include "adimensional_system_pcfm.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp_common/log.hpp" #include namespace specmicp { static constexpr int max_solver_restart = 3; // max number of time the solver is // restarted in case of failure // solve_equilibrium function // ========================== AdimensionalSystemSolution solve_equilibrium( std::shared_ptr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data, constraints, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } AdimensionalSystemSolution solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data_ptr, constraints, previous_solution, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } // The solver // ========== // constructor // ----------- AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution( Vector& x ) { set_true_variable_vector(x); return m_system->get_solution(x, m_var); } AdimensionalSystemSolution AdimensionalSystemSolver::unsafe_get_raw_solution( Vector& x ) { set_true_variable_vector(x); return m_system->unsafe_get_solution(x, m_var); } AdimensionalSystemSolution AdimensionalSystemSolver::merge_raw_solution_with_immobile_species( Vector& x, const AdimensionalSystemSolution& sol_immobile ) { set_true_variable_vector(x); AdimensionalSystemSolution sol = m_system->unsafe_get_solution(x, m_var); specmicp_assert(sol.main_variables.size() == sol_immobile.main_variables.size()); specmicp_assert(sol.sorbed_molalities.size() == sol_immobile.sorbed_molalities.size()); specmicp_assert(sol.log_gamma.size() == sol.log_gamma.size()); sol.sorbed_molalities = sol_immobile.sorbed_molalities; if (m_data->nb_mineral() > 0) { auto m_start = m_system->dof_mineral(0); sol.main_variables.segment(m_start, m_data->nb_mineral()) = sol_immobile.main_variables.segment(m_start, m_data->nb_mineral()); } if (m_data->nb_ssites() > 0) { auto s_start = m_system->dof_surface(0); sol.main_variables.segment(s_start, m_data->nb_ssites()) = sol_immobile.main_variables.segment(s_start, m_data->nb_ssites()); } return sol; } // Solving the system // ------------------ micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init) { // reset the units // - - - - - - - - m_system->set_units(get_options().units_set); // initialize the system // - - - - - - - - - - - if (init) { m_system->reasonable_starting_guess(x); if (get_options().use_pcfm) { run_pcfm(x); } } else if (get_options().force_pcfm) { run_pcfm(x); } // Solve the system // - - - - - - - - micpsolver::MiCPPerformance perf = solve_system(x); // If failed : try again // - - - - - - - - - - - int cnt = 0; while (perf.return_code < micpsolver::MiCPSolverReturnCode::Success and get_options().allow_restart and cnt < max_solver_restart ) { WARNING << "Failed to solve the system ! Return code :" << (int) perf.return_code << ". We shake it up and start again"; // Try safer options const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor; const bool save_scaling = get_options().solver_options.use_scaling; get_options().solver_options.use_scaling = true; if (save_penalization_factor == 1.0) get_options().solver_options.penalization_factor = 0.8; // Re-initialize the system set_return_vector(x); m_system->reasonable_restarting_guess(x); if (get_options().use_pcfm or get_options().force_pcfm) run_pcfm(x); // Solve the system micpsolver::MiCPPerformance perf2 = solve_system(x); // Restore options get_options().solver_options.penalization_factor = save_penalization_factor; get_options().solver_options.use_scaling = save_scaling; // Record the work performed perf += perf2; ++cnt; } // If solution is good // - - - - - - - - - - if (perf.return_code >= micpsolver::MiCPSolverReturnCode::Success) { set_return_vector(x); } return perf; } micpsolver::MiCPPerformance AdimensionalSystemSolver::solve_system(Vector& x) { set_true_variable_vector(x); if (not m_system->system_is_box_constrained()) { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_perfs(); } else { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_perfs(); } } // Variables management // --------------------- void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x) { // This function sets the true variable for the specmicp system // the order is important // this is synchronised with the specmicp system, but needs to be set here const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i {0}; if (m_system->is_active_component(0)) { m_var(new_i) = x(m_system->dof_water()); ++new_i; } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(), i ); if (it != non_active_component.cend()) continue; scalar_t value = x(m_system->dof_component(i)); if (value == -HUGE_VAL) // check for previously undefined value { value = m_system->get_options().new_component_concentration; } m_var(new_i) = value; ++new_i; } if (m_system->surface_model() != SurfaceEquationType::NoEquation) { for (index_t q: m_data->range_ssites()) { if (m_system->ideq_surf(q) != no_equation) { m_var(new_i) = x(m_system->dof_surface(q)); ++new_i; } } if (m_system->ideq_surf_pot() != no_equation) { m_var(new_i) = x(m_system->dof_surface_potential()); ++new_i; } } if (m_system->ideq_electron() != no_equation) { m_var(new_i) = x(m_system->dof_electron()); ++new_i; } - for (index_t m: m_data->range_mineral()) + if (m_system->get_options().solve_solid_solutions) { - if (m_system->is_active_mineral(m)) + for (auto d: m_data->range_solid_solutions()) { - m_var(new_i) = x(m_system->dof_mineral(m)); - ++new_i; + for (auto it=m_data->members_solid_solution(d);it;++it) + { + auto m=it.index(); + m_var(new_i) = x(m_system->dof_mineral(m)); + ++new_i; + } + } + for (auto m: m_data->range_mineral()) + { + if (m_system->is_active_mineral(m) and not m_data->is_mineral_in_solid_solution(m)) + { + m_var(new_i) = x(m_system->dof_mineral(m)); + ++new_i; + } + } + } else { + for (index_t m: m_data->range_mineral()) + { + if (m_system->is_active_mineral(m)) + { + m_var(new_i) = x(m_system->dof_mineral(m)); + ++new_i; + } } } m_var.conservativeResize(new_i); specmicp_assert(new_i == (unsigned) m_system->total_variables()); } void AdimensionalSystemSolver::set_return_vector(Vector& x) { // This function sets the variable from the specmicp "raw" solution // the order is important // this is synchronised with the specmicp system, but needs to be set here const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i = 0; if (m_system->is_active_component(0)) { x(m_system->dof_water()) = m_var(new_i); ++new_i; } else { x(m_system->dof_water()) = m_system->volume_fraction_water(x); } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(),i); if (it != non_active_component.cend()) { x(m_system->dof_component(i)) = -HUGE_VAL; continue; } x(m_system->dof_component(i)) = m_var(new_i) ; ++new_i; } for (index_t q: m_data->range_ssites()) { if (m_system->ideq_surf(q) != no_equation) { x(m_system->dof_surface(q)) = m_var(new_i); ++new_i; } else { x(m_system->dof_surface(q)) = -HUGE_VAL; } } if (m_system->ideq_surf_pot() != no_equation) { x(m_system->dof_surface_potential()) = m_var(new_i); ++new_i; } else { x(m_system->dof_surface_potential()) = -HUGE_VAL; } if (m_system->ideq_electron() != no_equation) { x(m_system->dof_electron()) = m_var(new_i); ++new_i; } else x(m_system->dof_electron()) = -HUGE_VAL; - for (index_t m: m_data->range_mineral()) + if (m_system->get_options().solve_solid_solutions) { - if (m_system->is_active_mineral(m)) + for (auto d: m_data->range_solid_solutions()) { - x(m_system->dof_mineral(m)) =m_var(new_i); - ++new_i; + for (auto it=m_data->members_solid_solution(d);it;++it) + { + auto m=it.index(); + x(m_system->dof_mineral(m)) = m_var(new_i); + ++new_i; + } } - else + for (auto m: m_data->range_mineral()) { - x(m_system->dof_mineral(m)) = 0.0; + if (m_system->is_active_mineral(m) and + (not m_data->is_mineral_in_solid_solution(m))) + { + x(m_system->dof_mineral(m)) = m_var(new_i); + ++new_i; + } + } + } else { + for (index_t m: m_data->range_mineral()) + { + if (m_system->is_active_mineral(m)) + { + x(m_system->dof_mineral(m)) = m_var(new_i); + ++new_i; + } + else + { + x(m_system->dof_mineral(m)) = 0.0; + } } } } // PCFM // ---- void AdimensionalSystemSolver::run_pcfm(Vector &x) { DEBUG << "Start PCFM initialization."; // we set up the true variable set_true_variable_vector(x); // The residual is computed to have some point of comparison Vector residuals(m_system->total_variables()); residuals.setZero(); m_system->get_residuals(m_var, residuals); const scalar_t res_0 = residuals.norm(); // the pcfm iterations are executed AdimensionalSystemPCFM pcfm_solver(m_system); PCFMReturnCode retcode = pcfm_solver.solve(m_var); // Check the answer if (retcode < PCFMReturnCode::Success) { // small prograss is still good enough m_system->get_residuals(m_var, residuals); const scalar_t final_residual = residuals.norm(); DEBUG << "Final pcfm residuals : " << final_residual << " set_secondary_variables(m_var); } } // Initialisation of variables // --------------------------- void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, std::unordered_map log_molalities, std::unordered_map volume_fraction_minerals, scalar_t log_free_sorption_site_concentration ) { m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; for (auto pair: log_molalities) { index_t idc = m_data->get_id_component(pair.first); if (idc == no_species or idc == m_data->electron_index() or idc == m_data->water_index()) { throw std::invalid_argument("This is not an aqueous component : "+pair.first); } if (pair.second > 0) { WARNING << "Initial molality for : " << pair.first << "is bigger than 1 molal."; } x(m_system->dof_component(idc)) = pair.second; } for (auto pair: volume_fraction_minerals) { index_t idm = m_data->get_id_mineral(pair.first); if (idm == no_species ) { throw std::invalid_argument("This is not a mineral at equilibrium : "+pair.first); } if (pair.second < 0 or pair.second > 1) { WARNING << "Initial volume fraction for : " << pair.first << "is not between 0 and 1."; } x(m_system->dof_mineral(idm)) = pair.second; } if (log_free_sorption_site_concentration != 0.0) { for (index_t q: m_data->range_ssites()) { x(m_system->dof_surface(q)) = log_free_sorption_site_concentration; } } } void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities ) { m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; if (log_molalities > 0) { WARNING << "Initial molality for : " << log_molalities << "is bigger than 1 molal."; } x.segment(1, m_data->nb_component()-1).setConstant(log_molalities); } void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities, scalar_t log_free_sorption_site_concentration ) { initialize_variables(x, volume_fraction_water, log_molalities); if (log_free_sorption_site_concentration != 0.0) { for (index_t q: m_data->range_ssites()) { x(m_system->dof_surface(q)) = log_free_sorption_site_concentration; } } } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solver_structs.cpp b/src/specmicp/adimensional/adimensional_system_solver_structs.cpp index 05356de..247aff4 100644 --- a/src/specmicp/adimensional/adimensional_system_solver_structs.cpp +++ b/src/specmicp/adimensional/adimensional_system_solver_structs.cpp @@ -1,65 +1,66 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "adimensional_system_solver_structs.hpp" #include "config_default_options_solver.h" namespace specmicp { AdimensionalSystemOptions::AdimensionalSystemOptions(): non_ideality(SPECMICP_DEFAULT_ENABLE_NONIDEAL), non_ideality_max_iter(SPECMICP_DEFAULT_NONIDEAL_MAX_ITER), scaling_electron(SPECMICP_DEFAULT_SCALING_ELECTRON), non_ideality_tolerance(SPECMICP_DEFAULT_NONIDEAL_TOL), under_relaxation_factor(SPECMICP_DEFAULT_UNDERRELAX_FACTOR), restart_concentration(SPECMICP_DEFAULT_RESTART_CONC), restart_water_volume_fraction(SPECMICP_DEFAULT_RESTART_WATER), new_component_concentration(SPECMICP_DEFAULT_NEW_COMPONENT_CONC), start_non_ideality_computation(SPECMICP_DEFAULT_TRHSOLD_START_NONIDEAL), - cutoff_total_concentration(SPECMICP_DEFAULT_CUTOFF_TOT_CONC) + cutoff_total_concentration(SPECMICP_DEFAULT_CUTOFF_TOT_CONC), + solve_solid_solutions(false) {} AdimensionalSystemSolverOptions::AdimensionalSystemSolverOptions(): allow_restart(SPECMICP_DEFAULT_ALLOW_RESTART), use_pcfm(SPECMICP_DEFAULT_USE_PCFM), force_pcfm(SPECMICP_DEFAULT_FORCE_PCFM) { units_set.length = SPECMICP_DEFAULT_LENGTH_UNIT; // disable the descent condition check solver_options.maxstep = SPECMICP_DEFAULT_MAX_STEP; solver_options.factor_descent_condition = SPECMICP_DEFAULT_FACTOR_DESC_COND; } } //end namespace spemicp diff --git a/src/specmicp/adimensional/adimensional_system_structs.hpp b/src/specmicp/adimensional/adimensional_system_structs.hpp index 3231c76..2b53691 100644 --- a/src/specmicp/adimensional/adimensional_system_structs.hpp +++ b/src/specmicp/adimensional/adimensional_system_structs.hpp @@ -1,478 +1,480 @@ /* ============================================================================= Copyright (c) 2014-2017 F. Georget Princeton University Copyright (c) 2017-2018 F. Georget 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 //! \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; + //! Solve solid solutions + bool solve_solid_solutions; 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 FixedSaturationIndex //!< Fixed Saturation Index }; //! \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) {} }; //! \struct FixedSaturationIndex //! \brief Struct to contain information needed to solve a fix saturation index struct FixedSaturationIndexConstraint { index_t id_component; //!< Index of the component index_t id_mineral; //!< Index of the mineral scalar_t saturation_index; //!< Value of the saturation index FixedSaturationIndexConstraint(index_t mineral, index_t component, scalar_t SI): id_component(component), id_mineral(mineral), saturation_index(SI) {} }; //! \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; //! \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}; //! Fixed fugacity constraints std::vector fixed_fugacity_cs {}; //! Fixed activity constraints std::vector fixed_activity_cs {}; //! Fixed molality constraints std::vector fixed_molality_cs {}; //! Fixed Saturation Index constraints std::vector fixed_SI_cs {}; //! Additional constraints for solid phases std::vector 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 Add a fixed saturation index for a mineral //! //! \param mineral id of the solid phase //! \param component id of the component //! \param saturation_index saturation index of the solid phase void add_fixed_saturation_index(index_t mineral, index_t component, scalar_t saturation_index) { fixed_SI_cs.emplace_back( FixedSaturationIndexConstraint(mineral, component, saturation_index) ); } //! \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/tests/specmicp/CMakeLists.txt b/tests/specmicp/CMakeLists.txt index e99a488..197cd4b 100644 --- a/tests/specmicp/CMakeLists.txt +++ b/tests/specmicp/CMakeLists.txt @@ -1,86 +1,88 @@ # Specmicp : Adim system # ---------------------- set(test_specmicp_adim_files adim/test_specmicp_adim.cpp adim/adimensional_system.cpp adim/adimensional_system_boxvi.cpp adim/adimensional_system_carboalu.cpp adim/adimensional_system_conditions.cpp adim/adimensional_system_configuration.cpp adim/adimensional_system_extractor.cpp adim/adimensional_system_oxydoreduc.cpp adim/adimensional_system_pcfm.cpp adim/adimensional_system_problem_solver.cpp + adim/adimensional_system_solidsolution.cpp adim/adimensional_system_solver.cpp adim/adimensional_system_sorption.cpp adim/adimensional_system_thermocarbo.cpp adim/adimensional_system_water_pressure.cpp adim/reactant_box.cpp ) if (HDF5_FOUND) list(APPEND test_specmicp_adim_files adim/io_hdf5_adimensional.cpp ) include_directories(HDF5_INCLUDE_DIRS) set_source_files_properties( adim/io_hdf5_adimensional.cpp PROPERTIES COMPILE_DEFINITIONS HDF5_DEFINITIONS ) endif() set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\") +set(TEST_CEMDATA18_PATH \"../../data/cemdata18.yaml\") set_source_files_properties(${test_specmicp_adim_files} PROPERTIES COMPILE_DEFINITIONS - "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}" + "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH};TEST_CEMDATA18_PATH=${TEST_CEMDATA18_PATH}" ) set(TEST_FEOH_PATH \"../../data/feoh_sorption.yaml\") set_property( SOURCE adim/adimensional_system_sorption.cpp APPEND PROPERTY COMPILE_DEFINITIONS TEST_FEOH_PATH=${TEST_FEOH_PATH} ) add_catch_test(NAME specmicp_adim SOURCES ${test_specmicp_adim_files} LINK_LIBRARIES ${SPECMICP_STATIC_LIBS} ) add_custom_target(adim_test_conf_file SOURCES adim/specmicp_conf_test.yaml) file(COPY adim/specmicp_conf_test.yaml DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) # Specmicp : Adim Kinetics # ------------------------ set(ADIMKINETICS_TEST_DIR ${PROJECT_TEST_DIR}/specmicp/adim_kinetics) set(test_specmicp_adim_kinetics_files ${ADIMKINETICS_TEST_DIR}/test_specmicp_adim_kinetics.cpp ${ADIMKINETICS_TEST_DIR}/kinetics_variables.cpp ${ADIMKINETICS_TEST_DIR}/set_mineral_kinetic.cpp ) set_source_files_properties(${test_specmicp_adim_kinetics_files} PROPERTIES COMPILE_DEFINITIONS "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}" ) add_catch_test(NAME specmicp_adim_kinetics SOURCES ${test_specmicp_adim_kinetics_files} LINK_LIBRARIES ${SPECMICP_STATIC_LIBS} ) # Specmicp : Binary # ----------------- add_subdirectory( binary ) diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index 786f010..9cf2808 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,755 +1,785 @@ #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 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 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 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 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 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 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 CEMDATA18") { + + std::cout << "\n Solving with reactant box and smart solver \n ------------------------ \n"; + + specmicp::database::Database thedatabase(TEST_CEMDATA18_PATH); + std::map swapping ({ + {"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(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 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 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 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 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 - fixed saturation index") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(1, "L/L"); reactant_box.add_fixed_saturation_index("Calcite", "HCO3[-]", 0.0); reactant_box.add_fixed_activity_component("H[+]", 1e-4); reactant_box.set_charge_keeper("Ca[2+]"); 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.pH() == Approx(4.0)); CHECK(extr.saturation_index(raw_data->get_id_mineral("Calcite")) == Approx(0.0).margin(1e-6)); } 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 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 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 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)); } } diff --git a/tests/specmicp/adim/adimensional_system_solidsolution.cpp b/tests/specmicp/adim/adimensional_system_solidsolution.cpp new file mode 100644 index 0000000..7a14964 --- /dev/null +++ b/tests/specmicp/adim/adimensional_system_solidsolution.cpp @@ -0,0 +1,332 @@ +#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 "specmicp_common/log.hpp" +#include + +static std::string mineral_appendix_database = R"plop( +[ +{ + "label": "Halite", + "composition": "Na[+], Cl[-]", + "log_k": -50, + "molar_volume": 50.0, + "flag_kinetic": 1 +} +] +)plop"; + + +static std::string ssol_appendix_database = R"plop( +[ +{ + "label": "AFmss", + "members": ["Monocarboaluminate", "Friedels_salt"] +} +] +)plop"; + +static std::string ssol_appendix_database18 = R"plop( +[ +{ + "label": "AFmss", + "members": ["Monocarbonate", "Friedels_salt"] +} +] +)plop"; + +static std::string ssol_CSHQ_appendix_database = R"plop( +[ +{ + "label": "CSHQss", + "members": ["CSHQ-JenD", "CSHQ-JenH"] #, "CSHQ-TobD", "CSHQ-TobH"] +} +] +)plop"; + + +TEST_CASE("Solid solution") +{ + specmicp::logger::ErrFile::stream() = &std::cerr; + specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; + + SECTION("Solid solution") + { + + std::cout << " -------------------------- \n" + << " Solid Solution (cemdata14)\n" + << " -------------------------- \n"; + specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); + thedatabase.swap_components({{"H[+]", "HO[-]"},{"Al[3+]", "Al(OH)4[-]"}}); + thedatabase.add_solid_phases(mineral_appendix_database, false); + thedatabase.add_solid_solutions(ssol_appendix_database, false); + 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, units_set); + + reactant_box.set_solution(15, "mL/L"); + reactant_box.add_solid_phase("Monocarboaluminate", 2.5, "g/L"); + reactant_box.add_solid_phase("Halite", 1.0, "g/L"); + reactant_box.add_aqueous_species("NaOH", 0.0, "mol/kg"); + + + auto constraints = reactant_box.get_constraints(true); + + specmicp::AdimensionalSystemSolverOptions opts; + opts.units_set = units_set; + opts.solver_options.maxstep = 10.0; + opts.solver_options.maxiter_maxstep = 500; + opts.solver_options.max_iter = 500; + opts.solver_options.use_crashing = false; + opts.solver_options.use_scaling = true; + opts.solver_options.set_tolerance(1e-8); + + opts.system_options.solve_solid_solutions = false; + + + + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, opts); + + + specmicp::Vector x; + solver.initialize_variables(x, 0.3, {{"HO[-]", -2}, {"Cl[-]", -1},}, {{"Monocarboaluminate", 0.1},{"Friedels_salt", 0.1}}); + auto perf = solver.solve(x, false); + + REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success); + + + specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + + specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set); + + std::cout << "pH : " << extr.pH() << "\n" + << "Hc : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Hemicarboaluminate")) + << " - " + << "Mc : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Monocarboaluminate")) + << " - " + << "Fs : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Friedels_salt")) + << "\n" + << "Hc : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Hemicarboaluminate"))*1e3 + << " - " + << "Mc : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Monocarboaluminate"))*1e3 + << " - " + << "Fs : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Friedels_salt"))*1e3 + << "\n" + << "CC : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Calcite"))*1e3 + << "\n" + << "CH : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Portlandite"))*1e3 + << "\n" + << "Cl[-] : " << extr.molality_component(raw_data->get_id_component("Cl[-]")) + << "\n" + << "Ca[2+] : " << extr.molality_component(raw_data->get_id_component("Ca[2+]")) + << "\n"; + + opts.system_options.solve_solid_solutions = true; + opts.solver_options.maxstep = 1.0; + + specmicp::AdimensionalSystemSolver solver2(raw_data, constraints, solution, opts); + specmicp::Vector x2 = solution.main_variables; + //x2(raw_data->nb_component()+raw_data->get_id_mineral("Hemicarboaluminate")) = 0.01; + auto perf2 = solver2.solve(x2, false); + REQUIRE(perf2.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success); + + specmicp::AdimensionalSystemSolution solution2 = solver2.get_raw_solution(x2);\ + specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, solver.get_options().units_set); + std::cout << "pH : " << extr2.pH() << "\n " + << "Hc : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Hemicarboaluminate")) + << " - " + << "Mc : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Monocarboaluminate")) + << " - " + << "Fs : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Friedels_salt")) + << "\n" + << "Hc : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Hemicarboaluminate"))*1e3 + << " - " + << "Mc : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Monocarboaluminate"))*1e3 + << " - " + << "Fs : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Friedels_salt"))*1e3 + << "\n" + << "CC : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Calcite"))*1e3 + << "\n" + << "CH : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Portlandite"))*1e3 + << "\n" + << "Cl[-] : " << extr2.molality_component(raw_data->get_id_component("Cl[-]")) + << "\n" + << "Ca[2+] : " << extr2.molality_component(raw_data->get_id_component("Ca[2+]")) + << "\n"; + } + + + SECTION("Solid solution") + { + std::cout << " -------------------------- \n" + << " Solid Solution (cemdata18)\n" + << " -------------------------- \n"; + specmicp::database::Database thedatabase(TEST_CEMDATA18_PATH); + thedatabase.swap_components({{"H[+]", "HO[-]"}}); + thedatabase.minerals_keep_only({"C4AH11", "Portlandite", "Monocarbonate", "Friedels_salt", "Hemicarbonate", "Calcite"}); + thedatabase.add_solid_phases(mineral_appendix_database, false); + thedatabase.add_solid_solutions(ssol_appendix_database18, false); + 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, units_set); + + reactant_box.set_solution(5, "mL/L"); + reactant_box.add_solid_phase("Monocarbonate", 2.5, "g/L"); + reactant_box.add_solid_phase("Halite", 1.0, "g/L"); + reactant_box.add_aqueous_species("NaOH", 0.0, "mol/kg"); + + + auto constraints = reactant_box.get_constraints(true); + + specmicp::AdimensionalSystemSolverOptions opts; + opts.units_set = units_set; + opts.solver_options.maxstep = 5.0; + opts.solver_options.maxiter_maxstep = 500; + opts.solver_options.max_iter = 500; + opts.solver_options.use_crashing = false; + opts.solver_options.use_scaling = true; + opts.solver_options.set_tolerance(1e-8); + + opts.system_options.solve_solid_solutions = false; + + + + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, opts); + + + specmicp::Vector x; + solver.initialize_variables(x, 0.5, {{"HO[-]", -2}, {"Cl[-]", -1}}, {{"Monocarbonate", 0.1},{"Friedels_salt", 0.1},{"C4AH11", 0.1}}); + auto perf = solver.solve(x, false); + + REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success); + + + specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + + specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set); + + std::cout << "pH : " << extr.pH() << "\n" + << "Hc : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Hemicarbonate")) + << " - " + << "Mc : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Monocarbonate")) + << " - " + << "Fs : " << extr.volume_fraction_mineral(raw_data->get_id_mineral("Friedels_salt")) + << "\n" + << "Hc : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Hemicarbonate"))*1e3 + << " - " + << "Mc : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Monocarbonate"))*1e3 + << " - " + << "Fs : " << extr.mass_concentration_mineral(raw_data->get_id_mineral("Friedels_salt"))*1e3 + << "\n" + << "Cl[-] : " << extr.molality_component(raw_data->get_id_component("Cl[-]")) + << "\n"; + + opts.system_options.solve_solid_solutions = true; + opts.solver_options.maxstep = 1.0; + solver = specmicp::AdimensionalSystemSolver(raw_data, constraints, solution, opts); + auto perf2 = solver.solve(x, false); + REQUIRE(perf2.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success); + + + specmicp::AdimensionalSystemSolution solution2 = solver.get_raw_solution(x); + + specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, solver.get_options().units_set); + std::cout << "pH : " << extr2.pH() << "\n " + << "Hc : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Hemicarbonate")) + << " - " + << "Mc : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Monocarbonate")) + << " - " + << "Fs : " << extr2.volume_fraction_mineral(raw_data->get_id_mineral("Friedels_salt")) + << "\n" + << "Hc : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Hemicarbonate"))*1e3 + << " - " + << "Mc : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Monocarbonate"))*1e3 + << " - " + << "Fs : " << extr2.mass_concentration_mineral(raw_data->get_id_mineral("Friedels_salt"))*1e3 + << "\n" + << "Cl[-] : " << extr2.molality_component(raw_data->get_id_component("Cl[-]")) + << "\n"; + } + +// SECTION("CSHQ Solid solution") +// { +// std::cout << "C-S-H-Q solid solution" << std::endl; + + +// specmicp::database::Database thedatabase(TEST_CEMDATA18_PATH); +// thedatabase.swap_components({{"H[+]", "HO[-]"}}); +// thedatabase.minerals_keep_only({"Portlandite", "CSHQ-JenD", "CSHQ-JenH"});//, "CSHQ-TobD", "CSHQ-TobH"}); +// thedatabase.add_solid_solutions(ssol_CSHQ_appendix_database, false); +// 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, units_set); + +// reactant_box.set_solution(100, "mL/L"); +// reactant_box.add_solid_phase("C3S", 200.0, "g/L"); +// //reactant_box.add_aqueous_species("NaOH", 0.1, "mol/kg"); + + +// auto constraints = reactant_box.get_constraints(true); + +// specmicp::AdimensionalSystemSolverOptions opts; +// opts.units_set = units_set; +// opts.solver_options.maxstep = 10.0; +// opts.solver_options.maxiter_maxstep = 500; +// opts.solver_options.max_iter = 500; +// opts.solver_options.use_crashing = false; +// opts.solver_options.use_scaling = true; +// opts.solver_options.set_tolerance(1e-8); + +// opts.system_options.solve_solid_solutions = false; + + + +// specmicp::AdimensionalSystemSolver solver(raw_data, constraints, opts); + + +// specmicp::Vector x; +// auto perf = solver.solve(x, true); + +// REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::Success); + + +// specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + +// specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set); + +// std::cout << extr.pH() << "\n " +// << "\n"; + +// opts.system_options.solve_solid_solutions = true; +// solver = specmicp::AdimensionalSystemSolver(raw_data, constraints, solution, opts); +// solver.solve(x, false); + +// specmicp::AdimensionalSystemSolution solution2 = solver.get_raw_solution(x); + +// specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, solver.get_options().units_set); +// std::cout << extr2.pH() << "\n " +// << "\n"; +// } +}