diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index 048f46d..10dffc2 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,2173 +1,2175 @@ /* ============================================================================= 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 // 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.get()), - m_equations(total_dofs(), ptrdata) + m_equations(total_dofs(), ptrdata), + m_constraints(constraints) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+ptrdata->nb_ssites()+1); - number_eq(constraints); + //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, ptrdata.get()), - m_equations(total_dofs(), ptrdata) + m_equations(total_dofs(), ptrdata), + m_constraints(constraints) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+ptrdata->nb_ssites()+1); - number_eq(constraints); + //number_eq(constraints); } // Secondary variables constructor // =============================== // No previous solution // -------------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( 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())), molar_fractions(Vector::Zero(data->nb_solid_solutions())), logactivity_solid_phases(Vector::Zero(data->nb_mineral())) {} // Previous solution // ----------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( 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), molar_fractions(Vector::Zero(data->nb_solid_solutions())), logactivity_solid_phases(Vector::Zero(data->nb_mineral())) { // 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_ssol(data->nb_solid_solutions(), true) {} // Equation numbering // ================== #define spc_tmp_HAS_SOLID_SOLUTION true #define spc_tmp_NO_SOLID_SOLUTION false template <> void AdimensionalSystem::number_eq_minerals( const AdimensionalSystemConstraints& constraints, index_t& neq ); template <> void AdimensionalSystem::number_eq_minerals( const AdimensionalSystemConstraints& constraints, index_t& neq ); template <> void AdimensionalSystem::number_eq_surfaces( const AdimensionalSystemConstraints& constraints, index_t& neq ); template <> void AdimensionalSystem::number_eq_surfaces( const AdimensionalSystemConstraints& constraints, index_t& neq ); // Note : this function also computes scaling factor that would be used in the computation // ------ void AdimensionalSystem::number_eq( - const AdimensionalSystemConstraints& constraints ) { + const AdimensionalSystemConstraints& constraints = m_constraints; + index_t neq = 0; // units compute_scaling_factors(); // Water // ===== if (constraints.water_equation != WaterEquationType::NoEquation) { number_eq_water(constraints, 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) { number_eq_surfaces(constraints, neq); } else if (constraints.surface_model.model_type == SurfaceEquationType::EDL) { number_eq_surfaces(constraints, neq); } else { m_equations.type_equation(dof_surface_potential()) = static_cast(SurfaceEquationType::NoEquation); } // Secondary species // ================= // Secondary aqueous species // ------------------------- bool solve_electron_equation = activate_secondary_species(constraints); // Gas // --- activate_gas(constraints); // 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 { activate_sorbed_species(constraints); } // Electron equation // ----------------- if (solve_electron_equation) { number_eq_electron(constraints, neq); } // Minerals // ======== if (constraints.immobile_disabled) { m_equations.nb_free_variables = neq; } else { if (get_options().solve_solid_solutions) { number_eq_minerals(constraints, neq); } else { number_eq_minerals(constraints, neq); } } m_equations.nb_tot_variables = neq; m_equations.nb_complementarity_variables = m_equations.nb_tot_variables - m_equations.nb_free_variables; } void AdimensionalSystem::number_eq_water( const AdimensionalSystemConstraints& constraints, index_t& neq ) { 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); } 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; } } } bool AdimensionalSystem::activate_secondary_species( const AdimensionalSystemConstraints& constraints ) { 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; } return solve_electron_equation; } void AdimensionalSystem::activate_gas( const AdimensionalSystemConstraints& constraints) { 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); } } void AdimensionalSystem::activate_sorbed_species( const AdimensionalSystemConstraints& constraints ) { // 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); } } } void AdimensionalSystem::number_eq_electron( const AdimensionalSystemConstraints& constraints, index_t& neq ) { 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; } } } } bool AdimensionalSystem::mineral_can_precipitate( index_t m, const AdimensionalSystemConstraints& constraints ) { 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; 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; } } } } // Not an equation if fixed saturation index for (const auto& it : constraints.fixed_SI_cs) { if (it.id_mineral == m) { can_precipitate = false; } } return can_precipitate; } template <> void AdimensionalSystem::number_eq_minerals( const AdimensionalSystemConstraints& constraints, index_t& neq ) { // above equations are 'free' (i.e. non constrained) // following equations are complementarity conditions m_equations.nb_free_variables = neq; for (index_t m: m_data->range_mineral()) { if (mineral_can_precipitate(m, constraints)) { m_equations.add_equation(dof_mineral(m), &neq); } } // Minerals => no precipitation // // /!\ must be set after equations ids are set if (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); } } } } template <> void AdimensionalSystem::number_eq_minerals( const AdimensionalSystemConstraints& constraints, index_t& neq ) { std::vector can_precipitate(m_data->nb_mineral(), true); for (index_t m: m_data->range_mineral()) { can_precipitate[m] = mineral_can_precipitate(m, constraints); } 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; } } // 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); } } // Minerals => no precipitation // // /!\ must be set after equations ids are set if (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); } } } } #undef spc_tmp_HAS_SOLID_SOLUTION #undef spc_tmp_NO_SOLID_SOLUTION template <> void AdimensionalSystem::number_eq_surfaces( const AdimensionalSystemConstraints& constraints, index_t& neq ) { 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; } template <> void AdimensionalSystem::number_eq_surfaces( const AdimensionalSystemConstraints& constraints, index_t& neq ) { 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); } } // 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) + 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) { 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(); 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(); 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 e9768c2..127e9a3 100644 --- a/src/specmicp/adimensional/adimensional_system.hpp +++ b/src/specmicp/adimensional/adimensional_system.hpp @@ -1,957 +1,963 @@ /* ============================================================================= 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 constraints Constraints to apply to the system (only reference is stored) //! \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 constraints Constraints to apply to the system (only reference is stored) //! \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 + //! + //! Only to be called before numbering the equations 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 // ################## +public: //! \brief Number the equations //! \param constraints the constraints to apply to the system //! \sa number_eq_aqueous_component - void number_eq(const AdimensionalSystemConstraints& constraints); + //void number_eq(const AdimensionalSystemConstraints& constraints); + void number_eq(); +private: //! \brief Number the equations for the water //! \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_water( const AdimensionalSystemConstraints& constraints, index_t& neq ); //! \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 ); //! \brief Check if secondary species need to be taken into account bool activate_secondary_species( const AdimensionalSystemConstraints& constraints ); //! \brief Check if gases need to be taken into account void activate_gas( const AdimensionalSystemConstraints& constraints ); //! \brief Check if sorbed species need to be taken into account void activate_sorbed_species( const AdimensionalSystemConstraints& constraints ); //! \brief Number the electron equation //! \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_electron( const AdimensionalSystemConstraints& constraints, index_t& neq ); //! \brief Number the equations for the minerals //! \param constraints the constraints to apply to the system //! \param[in,out] neq the number of equations in the system //! \sa number_eq template void number_eq_minerals( const AdimensionalSystemConstraints& constraints, index_t& neq ); //! \brief Check if a mineral can precipitate given the constraints bool mineral_can_precipitate( index_t m, const AdimensionalSystemConstraints& constraints ); //! \brief Number the equations for the surface model //! \param constraints the constraints to apply to the system //! \param[in,out] neq the number of equations in the system //! \sa number_eq template void number_eq_surfaces( 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(database::DataContainer* data); //! \brief Initialization with a 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; + const AdimensionalSystemConstraints& m_constraints; }; } // end namespace specmicp #endif // SPECMICP_ADIMENSIONALSYSTEM_HPP diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp index 702afc4..4b91de3 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.cpp +++ b/src/specmicp/adimensional/adimensional_system_solver.cpp @@ -1,594 +1,601 @@ /* ============================================================================= 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())) {} +AdimensionalSystemOptions AdimensionalSystemSolver::get_system_options() +{ + return m_system->get_options(); +} + 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_options(get_options().system_options); m_system->set_units(get_options().units_set); + m_system->number_eq(); // 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; } if (m_system->get_options().solve_solid_solutions) { for (auto d: m_data->range_solid_solutions()) { 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; if (m_system->get_options().solve_solid_solutions) { for (auto d: m_data->range_solid_solutions()) { 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; } } for (auto m: m_data->range_mineral()) { 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.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp index 636b0e0..f1db0f0 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.hpp +++ b/src/specmicp/adimensional/adimensional_system_solver.hpp @@ -1,238 +1,243 @@ /* ============================================================================= 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. * ============================================================================= */ #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #include "adimensional_system_solver_structs.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/options_handler.hpp" #include #include /*! \file adimensional_system_solver.hpp \brief Solve a reduced system \code{cpp} std::shared_ptr raw_data; AdimensionalSystemConstraints constraints; // ... AdimensionalSystemSolution previous_solution; // ... AdimensionalSystemSolverOptions options; // ... AdimensionalSystemSolver solver(raw_data, 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"); } AdimensionalSystemSolution solution = solver.get_raw_solution(variables); \endcode */ namespace specmicp { // forward declaration class AdimensionalSystem; struct AdimensionalSystemSolution; //! \brief Solve an adimensional system //! //! Take care of non-existing component in the system //! and restart the computation if necessary //! //! \ingroup specmicp_api class SPECMICP_DLL_PUBLIC AdimensionalSystemSolver: public OptionsHandler { public: //! \brief Pointer to the system using SystemPtr = std::shared_ptr; //! Default constructor //! //! Another constructor is necessary AdimensionalSystemSolver() {} //! \brief Initialise a solver from scratch //! //! \param data A raw database - //! \param constraints The system to solver + //! \param constraints The system to solve (only reference is stored) AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ); //! \brief Initialise a solver from scratch //! //! \param data A raw database - //! \param constraints The system to solver + //! \param constraints The system to solve (only reference is stored) //! \param options (optional) customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database - //! \param constraints The system to solver + //! \param constraints The system to solve (only reference is stored) //! \param previous_solution A solution of a similar system that will be used to initialize the system AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database - //! \param constraints The system to solver + //! \param constraints The system to solve (only reference is stored) //! \param previous_solution A solution of a similar system that will be used to initialize the system //! \param options customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); //! \brief solve the problem using initial guess x //! //! \param[in,out] x in -> initial guess, out -> solution //! \param init if true, the algorithm guess a starting point micpsolver::MiCPPerformance solve(Vector& x, bool init=false); //! \brief Return the system used for the computation SystemPtr get_system() {return m_system;} + //! \brief Return the system options + //! + //! /!\ be careful when changing this ! + AdimensionalSystemOptions get_system_options(); + //! \brief Return the solution in a manageable form //! //! \param x The solution (complete set of the main variables) AdimensionalSystemSolution get_raw_solution(Vector& x); //! \brief Return the solution in a manageable form //! //! \param x The solution (complete set of the main variables) AdimensionalSystemSolution unsafe_get_raw_solution(Vector& x); //! \brief Return the solution merge with solid phases of another //! //! This is useful when used with the "immobilized solid phases" AdimensionalSystemSolution merge_raw_solution_with_immobile_species( Vector& x, const AdimensionalSystemSolution& sol_immobile ); //! \brief Initialize the problem using the Positive continuous fraction method //! //! \sa specmicp::AdimensionalSystemPCFM void run_pcfm(Vector& x); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for chosen aqueous component //! \param volume_fraction_minerals volume fraction of the minerals //! \param log_free_sorption_site_concentration concentration of free sorption site void initialize_variables(Vector& x, scalar_t volume_fraction_water, std::unordered_map log_molalities, std::unordered_map volume_fraction_minerals = std::unordered_map(), scalar_t log_free_sorption_site_concentration = 0 ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components //! \param log_free_sorption_site_concentration concentration of free sorption site void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities, scalar_t log_free_sorption_site_concentration ); private: //! \brief set up the true variable vector //! //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_true_variable_vector(const Vector& x); //! \brief set up the true solution vector //! //! add zero components //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_return_vector(Vector& x); //! \brief solve the problem micpsolver::MiCPPerformance SPECMICP_DLL_LOCAL solve_system(Vector& x); RawDatabasePtr m_data; //! The raw database SystemPtr m_system; //! The system to solve Vector m_var; //! Copy of the solution vector (necessary in case of failing) }; //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints Constraints applied to the system //! \param options Options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints constraints applied to the system //! \param previous_solution a previous solution //! \param options options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); } // end namespace specmicp #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP diff --git a/tests/specmicp/adim/adimensional_system.cpp b/tests/specmicp/adim/adimensional_system.cpp index 89de992..93e3275 100644 --- a/tests/specmicp/adim/adimensional_system.cpp +++ b/tests/specmicp/adim/adimensional_system.cpp @@ -1,341 +1,348 @@ #include "catch.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp_database/database.hpp" #include "specmicp_database/appender.hpp" #include static std::string ssites_appendix_database = R"plop( [ { "label": "XS", } ] )plop"; static std::string sorbed_appendix_database = R"plop( [ { "label": "S1[2+]", "composition": "Ca[2+], XS", "log_k": -5 } ] )plop"; specmicp::RawDatabasePtr get_test_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_half_cell_reactions(); thedatabase.add_surface_sites(ssites_appendix_database); thedatabase.add_sorbed_species(sorbed_appendix_database); return thedatabase.get_database(); } specmicp::RawDatabasePtr get_test_oxydoreduc_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("Adimensional system", "[specmicp][MiCP program][adimensional]") { RawDatabasePtr thedatabase = get_test_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("initialization") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; AdimensionalSystemConstraints constraints(total_concentration); AdimensionalSystem system(thedatabase, constraints); + system.number_eq(); REQUIRE(system.ideq_w() != no_equation); REQUIRE(system.ideq_paq(id_oh) != no_equation); REQUIRE(system.ideq_paq(id_ca) != no_equation); REQUIRE(system.ideq_min(id_ch) != no_equation); Vector x = Vector::Zero(system.total_dofs()); x(id_h2o) = 1.0; x(id_oh) = std::log10(2e-3); x(id_ca) = -3; system.set_secondary_concentration(x); scalar_t molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); REQUIRE(molality_h/5e-12 == Approx(0.0).margin(1e-2)); system.compute_log_gamma(x); molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); system.hook_start_iteration(x, 1.0); molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); system.hook_start_iteration(x, 1.0); REQUIRE(molality_h == Approx(system.secondary_molality(0))); } SECTION("Residuals") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystem system(thedatabase, constraints); + system.number_eq(); Vector x = Vector::Zero(system.total_dofs()); x(system.ideq_w()) = 1.0; x(system.ideq_paq(id_oh)) = std::log10(2e-3); x(system.ideq_paq(id_ca)) = -3; specmicp::scalar_t res = system.residual_water_conservation(x); REQUIRE(std::isfinite(res)); REQUIRE(res == Approx((55.5 - system.density_water()/thedatabase->molar_mass_basis(id_h2o))/55.5)); } SECTION("Equation numbering") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.disable_conservation_water(); Vector sc(1); sc <<1.23456; constraints.set_equilibrium_surface_model(sc); constraints.total_concentrations(2) = 0; specmicp::AdimensionalSystem system(thedatabase, constraints); + system.number_eq(); CHECK(system.water_equation_type() == specmicp::WaterEquationType::NoEquation); CHECK(system.surface_model() == specmicp::SurfaceEquationType::Equilibrium); CHECK(system.surface_total_concentration(0) == 1.23456); CHECK((int) system.aqueous_component_equation_type(2) == (int) specmicp::AqueousComponentEquationType::NoEquation); constraints = specmicp::AdimensionalSystemConstraints(total_concentration); constraints.enable_conservation_water(); constraints.disable_surface_model(); constraints.add_fixed_activity_component(2, -2.0); - system = specmicp::AdimensionalSystem(thedatabase, constraints); + specmicp::AdimensionalSystem system2(thedatabase, constraints); + system2.number_eq(); - CHECK(system.water_equation_type() == specmicp::WaterEquationType::MassConservation); - CHECK(system.surface_model() == specmicp::SurfaceEquationType::NoEquation); - CHECK(system.aqueous_component_equation_type(2) == specmicp::AqueousComponentEquationType::FixedActivity); - CHECK(system.fixed_activity_bc(2) == -2.0); + CHECK(system2.water_equation_type() == specmicp::WaterEquationType::MassConservation); + CHECK(system2.surface_model() == specmicp::SurfaceEquationType::NoEquation); + CHECK(system2.aqueous_component_equation_type(2) == specmicp::AqueousComponentEquationType::FixedActivity); + CHECK(system2.fixed_activity_bc(2) == -2.0); } SECTION("Jacobian") { auto append = database::DataAppender(thedatabase); Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.enable_conservation_water(); Vector sc(1); sc <<1e-2; //constraints.set_equilibrium_surface_model(sc); constraints.set_EDL_surface_model(sc, 100); //constraints.add_fixed_activity_component(thedatabase->get_id_component("HO[-]"), -6); //constraints.set_charge_keeper(thedatabase->get_id_component("Ca[2+]")); - specmicp::AdimensionalSystem system (thedatabase, constraints); + specmicp::AdimensionalSystem system (thedatabase, constraints); + system.number_eq(); - Vector x = Vector::Zero(system.total_variables()); - x(system.ideq_w()) = 0.9; - x(system.ideq_paq(id_oh)) = std::log10(2e-3); - x(system.ideq_paq(id_ca)) = -3; - x(system.ideq_surf(0)) = -6; - x(system.ideq_surf_pot()) = 1e-2; + Vector x = Vector::Zero(system.total_variables()); + x(system.ideq_w()) = 0.9; + x(system.ideq_paq(id_oh)) = std::log10(2e-3); + x(system.ideq_paq(id_ca)) = -3; + x(system.ideq_surf(0)) = -6; + x(system.ideq_surf_pot()) = 1e-2; specmicp::Vector residual; system.get_residuals(x, residual); std::cout << residual << std::endl; specmicp::Matrix fd_jacobian; system.finite_difference_jacobian(x, fd_jacobian); specmicp::Matrix analytical_jacobian; system.analytical_jacobian(x, analytical_jacobian); for (index_t i: specmicp::range(thedatabase->nb_component() +thedatabase->nb_ssites() +thedatabase->nb_mineral()+1)) { std::cout << " } " << i << " -> " << system.ideq(i) << " } " << std::endl; } std::cout << thedatabase->nb_component() << " - " << thedatabase->nb_mineral() << " - " << thedatabase->nb_ssites() << " - " << std::endl; std::cout << analytical_jacobian - fd_jacobian << std::endl; specmicp::Matrix comp = (analytical_jacobian - fd_jacobian); CHECK(comp(system.ideq_paq(id_oh),system.ideq_surf_pot()) < 1e-2); CHECK(comp(system.ideq_surf_pot(),system.ideq_surf_pot()) < 1e-2); comp.col(system.ideq_surf_pot()).setZero(); std::cout << comp << std::endl; REQUIRE(comp.isMuchSmallerThan(1.0, 1e-3)); } SECTION("Jacobian - no surface") { auto append = database::DataAppender(thedatabase); Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.enable_conservation_water(); //Vector sc(1); //sc <<1.23456; //constraints.set_equilibrium_surface_model(sc); specmicp::AdimensionalSystem system (thedatabase, constraints); + system.number_eq(); Vector x = Vector::Zero(system.total_variables()); x(system.ideq_w()) = 1.0; x(system.ideq_paq(id_oh)) = std::log10(2e-3); x(system.ideq_paq(id_ca)) = -3; for (index_t i: specmicp::range(thedatabase->nb_component() +thedatabase->nb_ssites() +thedatabase->nb_mineral()+1)) { std::cout << " } " << i << " -> " << system.ideq(i) << std::endl; } specmicp::Vector residual; system.get_residuals(x, residual); std::cout << residual << std::endl; specmicp::Matrix analytical_jacobian; system.analytical_jacobian(x, analytical_jacobian); specmicp::Matrix fd_jacobian; system.finite_difference_jacobian(x, fd_jacobian); std::cout << thedatabase->nb_component() << " - " << thedatabase->nb_mineral() << " - " << thedatabase->nb_ssites() << " - " << std::endl; std::cout << analytical_jacobian - fd_jacobian << std::endl; REQUIRE((analytical_jacobian - fd_jacobian).isMuchSmallerThan(1.0, 1e-3)); } } TEST_CASE("Adimensional system - oxydoreduc", "[specmicp][MiCP program][adimensional]") { RawDatabasePtr thedatabase = get_test_oxydoreduc_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Jacobian - oxydoreduc") { // add a sorbed species, just to see if it works auto append = database::DataAppender(thedatabase); append.add_ssites("[{label: \"XS1\"}]"); append.add_sorbed(R"plop( - label: SorbedSpecies log_k: -1 composition: "2 HO[-], Ca[2+]" nb_sites_occupied: 1 )plop"); Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.enable_conservation_water(); specmicp::AdimensionalSystem system (thedatabase, constraints); + system.number_eq(); Vector x = Vector::Zero(system.total_variables()); x(system.ideq_w()) = 1.0; x(system.ideq_paq(id_oh)) = std::log10(2e-3); x(system.ideq_paq(id_ca)) = -3.0; x(system.ideq_electron()) = -2.0; for (index_t j: thedatabase->range_aqueous()){ std::cout << "\t" << thedatabase->get_label_aqueous(j); } std::cout << "\n"; std::cout << thedatabase->aqueous.get_nu_matrix() << std::endl; specmicp::Vector residual; system.get_residuals(x, residual); std::cout << residual << std::endl; specmicp::Matrix analytical_jacobian; system.analytical_jacobian(x, analytical_jacobian); specmicp::Matrix fd_jacobian; system.finite_difference_jacobian(x, fd_jacobian); std::cout << analytical_jacobian - fd_jacobian << std::endl; REQUIRE((analytical_jacobian - fd_jacobian).isMuchSmallerThan(1.0, 1e-3)); } } diff --git a/tests/specmicp/adim/adimensional_system_solver.cpp b/tests/specmicp/adim/adimensional_system_solver.cpp index e7c2b66..7db8d61 100644 --- a/tests/specmicp/adim/adimensional_system_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_solver.cpp @@ -1,251 +1,254 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/io/adimensional_system_solution_saver.hpp" #include "specmicp/io/adimensional_system_solution_reader.hpp" #include "specmicp_database/database.hpp" #include specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_half_cell_reactions(std::vector({"H2O", "HO[-]",})) ; return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("Solving adimensional system", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Solving simple case") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 3.0e4; total_concentration(id_oh) = 2.0e4; total_concentration(id_ca) = 1.0e4; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 100; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); + system->number_eq(); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.462142).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.820933).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volume in cm^3") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 10; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); units::UnitsSet unit_set; unit_set.length = units::LengthUnit::centimeter; system->set_units(unit_set); + system->number_eq(); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx( 0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volume in cm^3, mmol") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 1e3*0.03; total_concentration(id_oh) = 1e3*0.02; total_concentration(id_ca) = 1e3*0.01; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 10; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); units::UnitsSet unit_set; unit_set.length = units::LengthUnit::centimeter; unit_set.quantity = units::QuantityUnit::millimoles; system->set_units(unit_set); + system->number_eq(); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Automatic solver") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::Vector x; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints); solver.initialize_variables(x, 0.8, -2.0); //x(solver.dof_surface()) = -HUGE_VAL; solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; solver.get_options().solver_options.maxstep = 10.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.disable_descent_direction(); solver.get_options().solver_options.factor_gradient_search_direction = 100; auto perf = solver.solve(x); REQUIRE(perf.return_code >= micpsolver::MiCPSolverReturnCode::Success); auto unsafe_solution = solver.unsafe_get_raw_solution(x); auto solution = solver.get_raw_solution(x); AdimensionalSystemSolutionExtractor unsafe_extr( unsafe_solution, thedatabase, solver.get_options().units_set); AdimensionalSystemSolutionExtractor extr( solution, thedatabase, solver.get_options().units_set); // check that we have the good solution CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4)); CHECK(extr.log_molality_component(id_oh) == Approx(-1.46214).epsilon(1e-4)); CHECK(extr.log_molality_component(id_ca) == Approx(-1.82093).epsilon(1e-4)); //CHECK(extr.free_surface_concentration() == -HUGE_VAL); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4)); // check that unsafe solution is ok CHECK(extr.volume_fraction_water() == Approx(unsafe_extr.volume_fraction_water()).epsilon(1e-8)); CHECK(extr.log_molality_component(id_oh) == Approx(unsafe_extr.log_molality_component(id_oh)).epsilon(1e-8)); CHECK(extr.log_molality_component(id_ca) == Approx(unsafe_extr.log_molality_component(id_ca)).epsilon(1e-8)); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(unsafe_extr.volume_fraction_mineral(id_ch)).epsilon(1e-8)); CHECK(solution.log_gamma.norm() == Approx(unsafe_solution.log_gamma.norm()).epsilon(1e-8)); CHECK(solution.secondary_molalities.norm() == Approx(unsafe_solution.secondary_molalities.norm()).epsilon(1e-8)); database::Database the_db(thedatabase); the_db.save("test_db.yaml"); io::save_solution_yaml("test_solution.yaml", thedatabase, solution, "test_db.yaml"); RawDatabasePtr new_db; auto solution2 = io::parse_solution_yaml("test_solution.yaml", new_db); REQUIRE(new_db != nullptr); REQUIRE(solution.main_variables(0) == Approx(solution2.main_variables(0))); REQUIRE(solution.main_variables(id_oh) == Approx(solution2.main_variables(id_oh))); REQUIRE(solution.main_variables(id_ca) == Approx(solution2.main_variables(id_ca))); REQUIRE(extr.volume_fraction_mineral(id_ch) == Approx(solution2.main_variables(extr.dof_mineral(id_ch)))); REQUIRE(solution.log_gamma.norm() == Approx(solution2.log_gamma.norm())); CHECK(new_db->get_hash() == thedatabase->get_hash()); } }