diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index d4078c1..1e45852 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,1472 +1,1487 @@ /* ============================================================================= 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.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 // #define SPECMICP_DEBUG_EQUATION_FD_JACOBIAN namespace specmicp { //constexpr scalar_t log10f() const {return std::log(10.0);} //constexpr scalar_t log10 = log10f(); static const scalar_t log10 = std::log(10.0); // Constructor // =========== // No previous solution // -------------------- AdimensionalSystem::AdimensionalSystem( RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemOptions& options, const units::UnitsSet& units_set ): AdimemsionalSystemNumbering(ptrdata), OptionsHandler(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), m_second(ptrdata), m_equations(total_dofs(), ptrdata) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+1); number_eq(constraints); } // Previous solution // ----------------- AdimensionalSystem::AdimensionalSystem( RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemOptions& options, const units::UnitsSet& units_set ): AdimemsionalSystemNumbering(ptrdata), OptionsHandler(options), units::UnitBaseClass(units_set), m_inert_volume_fraction(constraints.inert_volume_fraction), m_second(previous_solution), m_equations(total_dofs(), ptrdata) { specmicp_assert(ptrdata->is_valid()); m_fixed_values.setZero(ptrdata->nb_component()+1); number_eq(constraints); } // Secondary variables constructor // =============================== // No previous solution // -------------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( const RawDatabasePtr& data ): secondary_molalities(Vector::Zero(data->nb_aqueous())), loggamma(Vector::Zero(data->nb_component()+data->nb_aqueous())), gas_fugacity(Vector::Zero(data->nb_gas())), gas_concentration(Vector::Zero(data->nb_gas())), sorbed_concentrations(Vector::Zero(data->nb_sorbed())) {} // Previous solution // ----------------- AdimensionalSystem::SecondaryVariables::SecondaryVariables( const AdimensionalSystemSolution& previous_solution ): secondary_molalities(previous_solution.secondary_molalities), loggamma(previous_solution.log_gamma), gas_fugacity(previous_solution.gas_fugacities), gas_concentration(Vector::Zero(previous_solution.gas_fugacities.rows())), sorbed_concentrations(previous_solution.sorbed_molalities) {} // IdEquations constructor // ======================= AdimensionalSystem::IdEquations::IdEquations( index_t nb_dofs, const RawDatabasePtr& data ): ideq(nb_dofs, no_equation), component_equation_type(data->nb_component()+1, no_equation), fixed_activity_species(data->nb_component()+1, no_species), active_aqueous(data->nb_aqueous(), false), active_gas(data->nb_gas(), false), active_sorbed(data->nb_sorbed()) {} // Equation numbering // ================== // Note : this function also computes scaling factor that would be used in the computation // ------ void AdimensionalSystem::number_eq( const AdimensionalSystemConstraints& constraints ) { index_t neq = 0; // units compute_scaling_factors(); // Water // ===== if (constraints.water_equation != WaterEquationType::NoEquation) { m_equations.type_equation(dof_water()) = static_cast(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; } } m_equations.add_equation(dof_water(), &neq); } // Aqueous components // ================== number_eq_aqueous_component(constraints, neq); // Surface model // ============= if (constraints.surface_model.model_type == SurfaceEquationType::Equilibrium) { // add the equation m_equations.add_equation(dof_surface(), &neq); m_equations.type_equation(dof_surface()) = static_cast(constraints.surface_model.model_type); // setup the total concentration m_fixed_values(dof_surface()) = constraints.surface_model.concentration; } // Secondary species // ================= // Secondary aqueous species // ------------------------- bool include_half_cell_reaction = (constraints.electron_constraint.equation_type != ElectronEquationType::NoEquation); bool solve_electron_equation {false}; for (auto j: m_data->range_aqueous()) { bool can_exist { true }; if ( include_half_cell_reaction or not m_data->is_half_cell_reaction(j)) for (const auto& k: m_equations.nonactive_component) { if (m_data->nu_aqueous(j, k) != 0.0) { can_exist = false; break; } } else { can_exist = false; } m_equations.set_aqueous_active(j, can_exist); if (can_exist and m_data->is_half_cell_reaction(j)) solve_electron_equation = true; //std::cout << m_data->labels_aqueous[j] << "can exist ? " << can_exist << " - logk : " << m_data->logk_aqueous(j) << std::endl; } // Gas // --- for (index_t k: m_data->range_gas()) { bool can_exist = true; for (const index_t& n: m_equations.nonactive_component) { if (m_data->nu_gas(k, n) != 0.0) { can_exist = false; break; } } m_equations.set_gas_active(k, can_exist); } // Sorbed species // -------------- for (index_t s: m_data->range_sorbed()) { // Check if the surface model is computed if (constraints.surface_model.model_type != SurfaceEquationType::Equilibrium) { m_equations.set_sorbed_active(s, false); continue; } // If so, check that all components of the sorbed species exist bool can_exist = true; for (const index_t& k: m_equations.nonactive_component) { if (m_data->nu_sorbed(s, k) != 0.0) { can_exist = false; break; } } m_equations.set_sorbed_active(s, can_exist); } // Electron equation // ----------------- if (solve_electron_equation) { m_equations.add_equation(dof_electron(), &neq); m_equations.type_equation(dof_electron()) = static_cast(constraints.electron_constraint.equation_type); if (constraints.electron_constraint.equation_type == ElectronEquationType::Equilibrium) { m_fixed_values(dof_electron()) = 0.0; } else if (constraints.electron_constraint.equation_type == ElectronEquationType::FixedpE) { m_fixed_values(dof_electron()) = constraints.electron_constraint.fixed_value; m_equations.related_species(dof_electron()) = constraints.electron_constraint.species; //assert(m_fixed_activity_species[dof_electron()] >= 0 // and m_fixed_activity_species[dof_electron()] < m_data->nb_aqueous()); //assert(m_data->is_half_cell_reaction(m_fixed_activity_species[dof_electron()])); } // scaling if (get_options().scaling_electron == 0.0) { for (index_t component : m_data->range_aqueous_component()) { if (aqueous_component_equation_type(component) == AqueousComponentEquationType::MassConservation) { get_options().scaling_electron = total_concentration_bc(component); break; } } } } // above equations are 'free' (i.e. non constrained) m_equations.nb_free_variables = neq; // following equations are complementarity conditions // Minerals // ======== for (index_t m: m_data->range_mineral()) { bool can_precipitate = true; // 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 } } if (can_precipitate) { m_equations.add_equation(dof_mineral(m), &neq); } } m_equations.nb_tot_variables = neq; m_equations.nb_complementarity_variables = m_equations.nb_tot_variables - m_equations.nb_free_variables; } 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; } // Finally number the equations for (index_t component: m_data->range_aqueous_component()) { // If no equation is assigned yet if (m_equations.type_equation(dof_component(component)) == static_cast(EqT::NoEquation)) { // Mass is conserved for this component //###FIXME: H[+], HO[-] const scalar_t& total_concentration = constraints.total_concentrations(dof_component(component)); if (std::abs(total_concentration) > get_options().cutoff_total_concentration) { m_equations.type_equation(dof_component(component)) = static_cast(EqT::MassConservation); m_fixed_values(dof_component(component)) = total_concentration; m_equations.add_equation(component, &neq); } else // add component to the nonactive component list { m_equations.add_non_active_component(component); } } // If equation is already assigned else { m_equations.add_equation(component, &neq); } } if (stdlog::ReportLevel() >= logger::Debug and m_equations.nonactive_component.size() > 0) { // if in debug mode list the non active components DEBUG << "Non active components :"; for (auto it: m_equations.nonactive_component) { DEBUG << " - " << it; } } } // Units // ================= void AdimensionalSystem::set_units(const units::UnitsSet& unit_set) { units::UnitBaseClass::set_units(unit_set); compute_scaling_factors(); } void AdimensionalSystem::compute_scaling_factors() { // /!\ scaling factors must be set at least once per call to this functions ! // the set should be = not *= because the scaling factors are typically set twice... // One at creation, the other before solving // (FIXME ?) (latter is necessary due to the facts that units can be changed) 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; } } // ================ // // // // Residuals // // // // ================ // 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 += log10*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::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 += log10*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 component) const { scalar_t sum = 0.0; for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; sum += log10*m_data->nb_sorption_sites(s)*m_data->nu_sorbed(s, component)*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 += log10*m_data->nu_gas(k, diff_component)*m_data->nu_gas(k, component)*gas_concentration(k); } return sum; } 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 (ideq_surf() != no_equation) 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 porosity = m_second.porosity; scalar_t sat_w = volume_fraction_water(x) / porosity; if (sat_w < 1) { const scalar_t pressure = m_equations.water_pressure_model(sat_w); res -= m_scaling_gas*(porosity - 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_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 (ideq_surf() != no_equation) aqconc += weigthed_sum_sorbed(component); res -= m_scaling_molality * conc_w * aqconc; res -= weigthed_sum_mineral(x, component); if (m_data->nb_gas() > 0) res -= weigthed_sum_gas(component); res /= total_concentration_bc(component); return res; } scalar_t AdimensionalSystem::residual_fixed_activity(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedActivity); scalar_t res = (fixed_activity_bc(component) - log_gamma_component(component) - log_component_molality(x, component) ); res /= fixed_activity_bc(component); return res; } scalar_t AdimensionalSystem::residual_fixed_molality(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedMolality); scalar_t res = (fixed_molality_bc(component) - log_component_molality(x, component) ); res /= fixed_molality_bc(component); return res; } scalar_t AdimensionalSystem::residual_fixed_fugacity(const Vector& x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedFugacity); index_t id_g = m_equations.fixed_activity_species[component]; scalar_t res = fixed_fugacity_bc(component) + m_data->logk_gas(id_g); for (index_t component: m_data->range_aqueous_component()) { if (m_data->nu_gas(id_g, component) == 0) continue; res -= m_data->nu_gas(id_g, component)*( log_gamma_component(component) + log_component_molality(x, component) ); } res /= fixed_fugacity_bc(component); return res; } scalar_t AdimensionalSystem::residual_mineral(const Vector& x, index_t m) const { specmicp_assert(ideq_min(m) != no_equation); scalar_t res = m_data->logk_mineral(m); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_mineral(m, i) != 0) { const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i); res -= m_data->nu_mineral(m, i)*log_activity_i; } } if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m)) res -= m_data->nu_mineral(m, m_data->electron_index())*log_activity_electron(x); return res; } scalar_t AdimensionalSystem::residual_charge_conservation(const Vector& x) const { scalar_t res = 0.0; for (index_t i: m_data->range_aqueous_component()) { if (m_data->charge_component(i) != 0 and ideq_paq(i) != no_equation) res += m_data->charge_component(i)*component_molality(x, i); } for (index_t j: m_data->range_aqueous()) { if (m_data->charge_aqueous(j) == 0 and not is_aqueous_active(j)) continue; res += m_data->charge_aqueous(j)*secondary_molality(j); } return m_scaling_molality*res; } scalar_t AdimensionalSystem::residual_surface(const Vector &x) const { specmicp_assert(ideq_surf() != no_equation); const scalar_t conc_w = density_water()*volume_fraction_water(x); scalar_t res = surface_total_concentration() - conc_w*free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; res -= m_scaling_molality*conc_w*m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } return res/surface_total_concentration(); } 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 (ideq_surf() != no_equation) 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 = 1 - sum_volume_fraction_minerals(x); set_volume_fraction_gas_phase(x); set_secondary_concentration(x); set_sorbed_concentrations(x); set_pressure_fugacity(x); // // water if (ideq_w() != no_equation) { switch (water_equation_type()) { case WaterEquationType::MassConservation: residual(ideq_w()) = residual_water_conservation(x); break; case WaterEquationType::SaturatedSystem: residual(ideq_w()) = residual_water_saturation(x); case WaterEquationType::NoEquation: break; } } // aqueous component for (index_t i: m_data->range_aqueous_component()) { switch (aqueous_component_equation_type(i)) { case AqueousComponentEquationType::NoEquation: break; case AqueousComponentEquationType::MassConservation: residual(ideq_paq(i)) = residual_component(x, i); break; case AqueousComponentEquationType::ChargeBalance: residual(ideq_paq(i)) = residual_charge_conservation(x); break; case AqueousComponentEquationType::FixedActivity: residual(ideq_paq(i)) = residual_fixed_activity(x, i); break; case AqueousComponentEquationType::FixedFugacity: residual(ideq_paq(i)) = residual_fixed_fugacity(x, i); break; case AqueousComponentEquationType::FixedMolality: residual(ideq_paq(i)) = residual_fixed_molality(x, i); break; } } // surface if (ideq_surf() != no_equation) residual(ideq_surf()) = residual_surface(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); // std::cout << residual << std::endl; } // ================ // // // // 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; 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 -= m_scaling_molality*diff_weigthed_sum_aqueous(k, 0); tmp -= m_scaling_molality*diff_weigthed_sum_sorbed(k, 0); // fixme gas tmp *= conc_w; tmp -= diff_weigthed_sum_gas(k, 0); jacobian(idw, ideq_paq(k)) = tmp/factor; } if (ideq_electron() != no_equation) { scalar_t tmp = 0.0; tmp -= m_scaling_molality*diff_weigthed_sum_aqueous(m_data->electron_index(), 0); tmp -= m_scaling_molality*diff_weigthed_sum_sorbed(m_data->electron_index(), 0); tmp*= conc_w; tmp -= diff_weigthed_sum_gas(m_data->electron_index(), 0); jacobian(idw, ideq_electron()) = tmp/factor; } 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; } } 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; } } } 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)*log10; tmp_iip -= diff_weigthed_sum_aqueous(k, i); tmp_iip -= 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)) = 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); tmp_iw -= weigthed_sum_aqueous(i); tmp_iw -= weigthed_sum_sorbed(i); tmp_iw *= m_scaling_molality*density_water(); jacobian(idp, ideq_w()) = tmp_iw/factor; } // Surface if (ideq_surf() != no_equation) { scalar_t tmp_s = -conc_w*m_scaling_molality*diff_surface_weigthed_sum_sorbed(i); jacobian(idp, ideq_surf()) = tmp_s/factor; } // Electron 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()) = 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); // 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; scalar_t tmp_value = m_data->nu_aqueous(j, k)*m_data->charge_aqueous(j); tmp_value *= secondary_molality(j)/component_molality(x, k); tmp_drdb += tmp_value; } jacobian(idp, idc) = m_scaling_molality*component_molality(x,k)*log10*tmp_drdb; } break; } // Fixed activity equation // ======================= case AqueousComponentEquationType::FixedActivity: { jacobian(idp, idp) = -1.0/fixed_activity_bc(i); break; } // Fixed fugacity equation // ======================= case AqueousComponentEquationType::FixedFugacity: { index_t id_g = m_equations.fixed_activity_species[i]; for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation or m_data->nu_gas(id_g, k) == 0.0) continue; jacobian(idp, ideq_paq(k)) = -m_data->nu_gas(id_g, k)/fixed_fugacity_bc(i); } break; } // end case // Fixed molality component // ======================== case AqueousComponentEquationType::FixedMolality: { jacobian(idp, idp) = -1.0/fixed_molality_bc(i); break; } } // end switch } } void AdimensionalSystem::jacobian_minerals(Vector& x, Matrix& jacobian) { for (index_t m: m_data->range_mineral()) { const index_t idm = ideq_min(m); if (idm == no_equation) continue; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) continue; jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i); } if (ideq_electron() != no_equation and m_data->is_mineral_half_cell_reaction(m)) jacobian(idm, ideq_electron()) = -m_data->nu_mineral(m, m_data->electron_index()); } } void AdimensionalSystem::jacobian_surface(Vector& x, Matrix& jacobian) { const index_t ids = ideq_surf(); const scalar_t factor = surface_total_concentration(); const scalar_t conc_w = density_water()*volume_fraction_water(x); specmicp_assert(ids != no_equation); scalar_t tmp_s = - free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; tmp_s -= m_data->nb_sorption_sites(s)* m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } jacobian(ids, ids) = m_scaling_molality*conc_w*log10 * tmp_s / factor; // water const index_t idw = ideq_w(); if (idw != no_equation) { const scalar_t rho_w = density_water(); scalar_t tmp_w = - free_sorption_site_concentration(x); for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) continue; tmp_w -= m_data->nb_sorption_sites(s)*sorbed_species_concentration(s); } jacobian(ids, idw) = m_scaling_molality*rho_w * tmp_w / factor; } // component for (index_t k: m_data->range_aqueous_component()) { const index_t idk = ideq_paq(k); if (idk == no_equation) continue; scalar_t tmp_k = - conc_w*diff_surface_weigthed_sum_sorbed(k); jacobian(ids, idk) = tmp_k/factor; } if (ideq_electron() != no_equation) { scalar_t tmp_e = - conc_w*diff_surface_weigthed_sum_sorbed(m_data->electron_index()); jacobian(ids, ideq_electron()) = tmp_e/factor; } // water } 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)) = 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) = 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 (ideq_surf() != no_equation) set_sorbed_concentrations(x); set_secondary_concentration(x); compute_log_gamma(x); set_volume_fraction_gas_phase(x); set_pressure_fugacity(x); } void AdimensionalSystem::set_volume_fraction_gas_phase(const Vector& x) { m_second.volume_fraction_gas = m_second.porosity - volume_fraction_water(x); } void AdimensionalSystem::set_pressure_fugacity(const Vector& x) { const auto rt = constants::gas_constant*temperature(); for (index_t k: m_data->range_gas()) { if (not is_active_gas(k)) continue; scalar_t logp = -m_data->logk_gas(k); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_gas(k, i) == 0.0) continue; const auto log_activity_i = log_component_molality(x, i) + log_gamma_component(i); logp += m_data->nu_gas(k, i) * log_activity_i; } if (ideq_electron() != no_equation and m_data->is_gas_half_cell_reaction(k)) logp += m_data->nu_gas(k, m_data->electron_index())*log_activity_electron(x); m_second.gas_fugacity(k) = pow10(logp); const scalar_t pressure = gas_fugacity(k)*gas_total_pressure(); const scalar_t concentration = m_scaling_gas*volume_fraction_gas_phase()*pressure/rt; m_second.gas_concentration(k) = concentration; } } void AdimensionalSystem::set_secondary_concentration(const Vector& x) { for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) { m_second.secondary_molalities(j) = 0.0; continue; } scalar_t logconc = -m_data->logk_aqueous(j) - log_gamma_secondary(j); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_aqueous(j, k) == 0) continue; const auto log_activity_k = log_component_molality(x, k) + log_gamma_component(k); logconc += m_data->nu_aqueous(j, k)*log_activity_k; } if (ideq_electron() != no_equation and m_data->is_half_cell_reaction(j)) logconc += m_data->nu_aqueous(j, m_data->electron_index())*log_activity_electron(x); m_second.secondary_molalities(j) = pow10(logconc); } } void AdimensionalSystem::set_sorbed_concentrations(const Vector& x) { for (index_t s: m_data->range_sorbed()) { if (not is_active_sorbed(s)) { m_second.sorbed_concentrations(s) = 0.0; continue; } scalar_t logconc = -m_data->logk_sorbed(s) + m_data->nb_sorption_sites(s)*(log_free_sorption_site_concentration(x)); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_sorbed(s, k) == 0.0) continue; const auto activity_k = log_component_molality(x, k) + log_gamma_component(k); logconc += m_data->nu_sorbed(s, k)*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); } } 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()); - for (index_t i: m_data->range_aqueous_component()) + + vector_debye_huckel_component(ion_str, sqrti); + vector_debye_huckel_aqueous(ion_str, sqrti); +} + +bool 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()) { - if (ideq_paq(i) == no_equation) - { - log_gamma_component(i) = 0.0; - continue; - } - log_gamma_component(i) = laws::extended_debye_huckel( - ionic_strength(), sqrti, - m_data->charge_component(i), - m_data->a_debye_component(i), - m_data->b_debye_component(i) - ); + 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; } - for (index_t j: m_data->range_aqueous()) +} + + +bool 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()) { - if (not is_aqueous_active(j)) - { - log_gamma_secondary(j) = 0.0; - continue; - } - log_gamma_secondary(j) = laws::extended_debye_huckel( - ionic_strength(), sqrti, - m_data->charge_aqueous(j), - m_data->a_debye_aqueous(j), - m_data->b_debye_aqueous(j) - ); + const scalar_t zisquare = std::pow(ionic_param(j, IoParam::charge_ind), 2); + log_gamma(offset+j) = - (constants::Adebye * zisquare * sqrt_ionic); + log_gamma(offset+j) /= ( 1 + ionic_param(j, IoParam::adebye_ind)*constants::Bdebye*sqrt_ionic); + log_gamma(offset+j) += ionic_param(j, IoParam::bdebye_ind)*ionic_strength; } } bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual) { if (not get_options().non_ideality) { set_secondary_variables(x); // we still need to compute secondary species ! // if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) // { // set_secondary_variables(x); // } return true; } not_in_linesearch = true; scalar_t previous_norm = m_second.loggamma.norm(); if (previous_norm == 0) previous_norm = 1.0; bool may_have_converged = false; if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) { m_equations.solve_pressure_model = true; // Use fixed point iterations for non-ideality for (int i=0; i 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::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 (ideq_surf() != no_equation) xtot(dof_surface()) = std::log10(0.8*surface_total_concentration()); else xtot(dof_surface()) = -HUGE_VAL; if (ideq_electron() != no_equation) xtot(dof_electron()) = -4; else xtot(dof_electron()) = -HUGE_VAL; xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero(); m_second = SecondaryVariables(m_data); } void AdimensionalSystem::reasonable_restarting_guess(Vector& xtot) { std::random_device rd; std::mt19937 gen(rd()); 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 (ideq_surf() != no_equation) xtot(dof_surface()) = std::log10(0.8*surface_total_concentration()); else xtot(dof_surface()) = -HUGE_VAL; if (ideq_electron() != no_equation) xtot(dof_electron()) = -4; else xtot(dof_electron()) = -HUGE_VAL; xtot.segment(offset_minerals(), m_data->nb_mineral()).setZero(); m_second = SecondaryVariables(m_data); } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp index bcb7980..380d1f6 100644 --- a/src/specmicp/adimensional/adimensional_system.hpp +++ b/src/specmicp/adimensional/adimensional_system.hpp @@ -1,794 +1,803 @@ /* ============================================================================= 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_ADIMENSIONALSYSTEM_HPP #define SPECMICP_ADIMENSIONALSYSTEM_HPP //! \file adimensional_system.hpp The MiCP program to solve speciation #include "adimensional_system_numbering.hpp" #include "adimensional_system_structs.hpp" #include "specmicp_common/types.hpp" #ifndef SPECMICP_DATABASE_HPP #include "specmicp_database/database.hpp" #endif #include "specmicp_common/micpsolver/micpprog.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/options_handler.hpp" //! \namespace specmicp main namespace for the SpecMiCP solver namespace specmicp { struct AdimensionalSystemSolution; //! \brief The equilibrium speciation problem //! //! Represent the equilibrium speciation problem as a Mixed Complementarity program //! //! The main variables are : //! - \f$phi_w\f$ the volume fraction of water (total saturation) //! - \f$b_i\f$ the molality of aqueous component //! - \f$phi_m\f$ the volume fraction of mineral //! - \f$s_f\f$ the concentration (molality) of free sorption sites //! //! The solid phase equilibrium is solved by using a complementarity formulation //! \f[ S^t_m \geq 0, \; - \mathrm{SI}_m \geq 0 , \; \mathrm{and} \; - S^t_m \mathrm{SI}_m = 0 \f] //! //! The secondary variables (molalities of secondary aqueous species, gas fugacities, ionic strength, ...) //! are solved using a fixed point iteration at the beginning of each iteration. //! //! This class should only be used through the AdimensionalSystemSolver. class SPECMICP_DLL_LOCAL AdimensionalSystem : public micpsolver::MiCPProg, public AdimemsionalSystemNumbering, public OptionsHandler, private units::UnitBaseClass { public: //! \brief Create a reduced system //! //! \param ptrdata Shared pointer to the thermodynamic database //! \param constraints Constraints to apply to the system //! \param options Numerical options, mainly related to the secondary variables //! \param units_set The set of units AdimensionalSystem(RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemOptions& options=AdimensionalSystemOptions(), const units::UnitsSet& units_set=units::UnitsSet() ); //! \brief Create a reduced system, initialize the system using a previous solution //! //! \param ptrdata Shared pointer to the thermodynamic database //! \param constraints Constraints to apply to the system //! \param previous_solution The previous solution to use for initialization //! \param options Numerical options, mainly related to the secondary variables //! \param units_set The set of units AdimensionalSystem( RawDatabasePtr& ptrdata, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemOptions& options=AdimensionalSystemOptions(), const units::UnitsSet& units_set=units::UnitsSet() ); // Variables // ========== //! \name Variables //! \brief How many ? //! //! @{ //! \brief Return the total number of variables index_t total_variables() const {return m_equations.nb_tot_variables;} //! \brief Return the number of 'free' variables (not subject to complementarity conditions) index_t nb_free_variables() const {return m_equations.nb_free_variables;} //! \brief Return the number of variables subjected to the complementarity conditions index_t nb_complementarity_variables() const {return m_equations.nb_complementarity_variables;} //! @} //! \brief Return a pointer to the database RawDatabasePtr get_database() {return m_data;} // The linear system // ================== //! \name Residuals and Jacobian //! \brief Compute the residuals and the jacobian //! //! @{ //! \brief Return the residuals //! //! \param x Vector of the main variables //! \param[out] residual Vector containing the residuals void get_residuals(const Vector& x, Vector& residual); //! \brief Return the jacobian //! //! \param x Vector of the main variables //! \param[out] jacobian Dense matrix representing the jacobian void get_jacobian(Vector& x, Matrix& jacobian); //! \brief Return the residual for the water conservation equation //! //! \param x Vector of the main variables scalar_t residual_water_conservation(const Vector& x) const; //! \brief Return the residual for the water saturation equation //! //! \param x Vector of the main variables scalar_t residual_water_saturation(const Vector& x) const; //! \brief Return the residual for the 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) const; //! \brief Equilibrium condition for the minerals //! \param x Vector of the main variables //! \param m Index of the mineral (in the database) scalar_t residual_mineral(const Vector& x, index_t m) const; //! \brief Return the residual of the charge conservation equation //! //! This equation may be used instead of a mass conservation equation to //! explicitely enforce the electroneutrality. //! The component corresponding to this equation is called the charge keeper //! \param x Vector of the main variables scalar_t residual_charge_conservation(const Vector& x) const; //! \brief Return the residual corresponding to a fixed fugacity equation //! \param x Vector of the main variables //! \param component Index of the corresponding component (in the database) scalar_t residual_fixed_fugacity(const Vector& x, index_t component) const; //! \brief Return the residual corresponding to a fixed activity equation //! //! If 'component ' is charged, then a charge keeper should be set. //! \param x Vector of the main variables //! \param component Index of the corresponding component (in the database) scalar_t residual_fixed_activity(const Vector& x, index_t component) const; //! \brief Return the residual corresponding to a fixed molality equation //! //! If 'component ' is charged, then a charge keeper should be set. //! \param x Vector of the main variables //! \param component Index of the corresponding component (in the database) scalar_t residual_fixed_molality(const Vector& x, index_t component) const; //! \brief Return the residual corresponding to the electron equation scalar_t residual_electron(const Vector &x) const; //! \brief Compute the jacobian analytically void analytical_jacobian(Vector& x, Matrix& jacobian); //! \brief Compute the jacobian using finite difference void finite_difference_jacobian(Vector& x, Matrix& jacobian); //! \brief Return a scale factor to avoid negative mass during Newton's iteration //! //! \param x Vector of the main variables //! \param update Update of the current iteration double max_lambda(const Vector &x, const Vector &update); //! @} // Getters // ####### // Equation id access // ================== //! \name Equation Id //! \brief The equation id is the row of the equation in the residuals/jacobian //! //! The degree of freedom getter function are inherited from specmicp::AdimensionalSystemNumbering //! @{ //! \brief Return the equation id //! \param i Index of the main variable (in the set of the main variables) index_t ideq(index_t i) const {return m_equations.ideq[i];} //! \brief Return the equation number for conservation of water index_t ideq_w() const {return m_equations.ideq[dof_water()];} //! \brief Return the equation number of component 'i' //! \param i Index of the component (in the database) index_t ideq_paq(index_t i) const { specmicp_assert(i < m_data->nb_component()); return m_equations.ideq[dof_component(i)]; } //! \brief Return the equation number of the electron equation index_t ideq_electron() const { return m_equations.ideq[dof_electron()]; } //! \brief Return the equation number of the free surface concentration index_t ideq_surf() const { return m_equations.ideq[dof_surface()]; } //! \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 { if (ideq_surf() == no_equation) return SurfaceEquationType::NoEquation; return SurfaceEquationType::Equilibrium; } //! \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 scalar_t density_water() 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 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) const { return pow10(log_free_sorption_site_concentration(x)); } //! \brief Return the log_10 of the free sorption site concentration scalar_t log_free_sorption_site_concentration(const Vector& x) const { specmicp_assert(ideq_surf() != no_equation); return x(ideq_surf()); } //! \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 activity of 'component' //! //! Only use this method if 'component' has a fixed activity constraint //! \param component Index of the aqueous component (in the database) scalar_t fixed_activity_bc(index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedActivity); return m_fixed_values(component); } //! \brief Return the fixed molality of 'component' //! //! Only use this method if 'component' has a fixed molality constraint //! \param component Index of the aqueous component (in the database) scalar_t fixed_molality_bc(index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedMolality); return m_fixed_values(component); } //! \brief Return the fixed fugacity value for 'component' scalar_t fixed_fugacity_bc(index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::FixedFugacity); return m_fixed_values(component); } //! \brief Return the total concentration for the electron scalar_t electron_total_concentration() const { return 0.0; } // Gas // --- //! \brief Return the fugacity of 'gas' //! //! \param gas Index of the gas (in the database) scalar_t gas_fugacity(index_t gas) const { return m_second.gas_fugacity(gas); } //! \brief Return the concentration of 'gas' in the system //! //! \param gas Index of the gas (in the database) scalar_t gas_concentration(index_t gas) const { return m_second.gas_concentration(gas); } //! \brief Return true if gas is active //! //! \param gas Index of the gas (in the database) bool is_active_gas(index_t gas) const { return m_equations.active_gas[gas]; } //! \brief Return the volume fraction (total saturation) of the gas phase scalar_t volume_fraction_gas_phase() const noexcept { return m_second.volume_fraction_gas; } // Sorbed species // -------------- //! \brief Return the surface total concentration const scalar_t& surface_total_concentration() const { return m_fixed_values(dof_surface()); } //! \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.16+25; } //! @} // Solution // ================= //! \brief Return the equilibrium state of the system, the Solution of the speciation problem //! //! Also checks that the secondary variables were correctly solved //! //! \param xtot the complete set of main variables //! \param x the reduced set of main variables, only the variables with an equation AdimensionalSystemSolution get_solution(Vector& xtot, const Vector& x); //! \brief Return the equilibrium state of the system, the Solution of the speciation problem //! //! This version does not check anything. //! //! \param xtot the complete set of main variables //! \param x the reduced set of main variables, only the variables with an equation AdimensionalSystemSolution unsafe_get_solution(Vector& xtot, const Vector& x); //! \name Secondary variables computation //! \brief Algorithms to compute the secondary variables //! //! @{ //! \brief Compute the ionic strength //! //! \param x Vector of the main variables void set_ionic_strength(const Vector& x); //! \brief Compute the activity coefficients //! //! \param x Vector of the main variables void compute_log_gamma(const Vector& x); //! \brief Compute the secondary aqueous species molalities //! //! \param x Vector of the main variables void set_secondary_concentration(const Vector& x); //! \brief Compute the secondary variables //! //! \param x Vector of the main variables void set_secondary_variables(const Vector& x); //! \brief Compute the gas phase volume fraction //! //! \param x Vector of the main variables void set_volume_fraction_gas_phase(const Vector& x); //! \brief Compute the fugacity for all the gas //! //! \param x Vector of the main variables void set_pressure_fugacity(const Vector& x); //! \brief Compute the sorbed species concentrations //! //! \param x Vector of the main variables void set_sorbed_concentrations(const Vector& x); //! \brief This function is called at the beginning of each iteration by the solver //! //! It does the fixed-point iteration to compute the secondary variables //! \param x Vector of the main variables //! \param norm_residual Norm of the current residuals bool hook_start_iteration(const Vector &x, scalar_t norm_residual); //! @} //! \brief A reasonable (well... maybe...) starting guess //! //! \param xtot Vector of the main variables //! \sa #reasonable_starting_guess void reasonable_starting_guess(Vector& xtot); //! \brief A reasonable (maybe...) restarting guess //! //! \param xtot Vector of the main variables //! \sa #reasonable_restarting_guess void reasonable_restarting_guess(Vector& xtot); //! \brief Set the units void set_units(const units::UnitsSet& unit_set); // Private Members and attributes // ############################## private: + bool vector_debye_huckel_component( + scalar_t ionic_strength, + scalar_t sqrt_ionic + ); + bool 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 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 the free surface concentration //! \sa weigthed_sum_sorbed scalar_t diff_surface_weigthed_sum_sorbed(index_t component) 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 te 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 electron equation to te Jacobian //! \param x Reduced vector of the main variables //! \param jacobian Dense matrix containing the jacobian void jacobian_electron(Vector& x, Matrix& jacobian); // Equation numbering // ################## //! \brief Number the equations //! \param constraints the constraints to apply to the system //! \sa number_eq_aqueous_component void number_eq(const AdimensionalSystemConstraints& constraints); //! \brief Number the equations for the aqueous components //! \param constraints the constraints to apply to the system //! \param[in,out] neq the number of equations in the system //! \sa number_eq void number_eq_aqueous_component( const AdimensionalSystemConstraints& constraints, index_t& neq ); // Setter // ###### // main variables // -------------- // they are set by the solver // Secondary variables // ------------------- //! \brief Return a reference to log_10(γ_i) where i is a component scalar_t& log_gamma_component(index_t component) { return m_second.loggamma(dof_component_gamma(component)); } //! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species scalar_t& log_gamma_secondary(index_t aqueous) { return m_second.loggamma(dof_aqueous_gamma(aqueous)); } //! \brief Return a reference to the ionic strength scalar_t& ionic_strength() { return m_second.ionic_strength; } // Attributes // ########## //! \struct SecondaryVariables //! \brief Contains information about the secondary variables struct SecondaryVariables { scalar_t ionic_strength {0.0}; //!< The ionic Strength scalar_t volume_fraction_gas {0.0}; //!< The gas saturation scalar_t porosity {0.0}; //!< The porosity Vector secondary_molalities; //!< The secondary molalities Vector loggamma; //!< The log of activity coefficients Vector gas_fugacity; //!< The gas fugacities Vector gas_concentration; //!< The gas concentrations Vector sorbed_concentrations; //!< The sorbed concentrations //! \brief Initialization without a previous solution SecondaryVariables(const RawDatabasePtr& data); //! \brief Initialization with a previous solution SecondaryVariables(const AdimensionalSystemSolution& previous_solution); }; //! \struct IdEquations //! \brief BookKeeper for the equations id and type struct IdEquations { bool use_water_pressure_model {false}; water_partial_pressure_f water_pressure_model {nullptr}; bool solve_pressure_model {false}; index_t nb_tot_variables {0}; index_t nb_free_variables {0}; index_t nb_complementarity_variables {0}; std::vector 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; //! \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; } }; bool not_in_linesearch {true}; scalar_t m_inert_volume_fraction; scalar_t m_scaling_molar_volume {1.0}; scalar_t m_scaling_gas {1.0}; scalar_t m_scaling_molality {1.0}; Vector m_fixed_values {}; SecondaryVariables m_second; IdEquations m_equations; }; } // end namespace specmicp #endif // SPECMICP_ADIMENSIONALSYSTEM_HPP diff --git a/src/specmicp_database/species/aqueous_list.hpp b/src/specmicp_database/species/aqueous_list.hpp index f3af948..e6c389f 100644 --- a/src/specmicp_database/species/aqueous_list.hpp +++ b/src/specmicp_database/species/aqueous_list.hpp @@ -1,175 +1,181 @@ /* ============================================================================= 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_DATABASE_AQUEOUSLIST_HPP #define SPECMICP_DATABASE_AQUEOUSLIST_HPP #include "specmicp_common/types.hpp" #ifndef SPECMICP_DATABASE_SPECIES_HPP #include "species.hpp" #endif #ifndef SPECMICP_DATABASE_IONICPARAMETERS_HPP #include "ionic_parameters.hpp" #endif //! \file aqueous_list.hpp //! \brief The list of aqueous species namespace specmicp { + +class AdimensionalSystem; + namespace database { //! \struct AqueousValues //! \brief Struct to initialize an aqueous species in the aqueous species list //! //! \ingroup database_species struct AqueousValues { std::string label; scalar_t logk; IonicModelValues ionic_values; }; //! \class AqueousList //! \brief The aqueous species //! //! This is the list of aqueous species in the system //! //! \ingroup database_species class AqueousList: public ReactiveSpeciesList { public: //! Initialize an empty list AqueousList() {} //! Initialize a list of 'siz' aqueous species with 'nb_components' components AqueousList(index_t size, index_t nb_components): ReactiveSpeciesList(size, nb_components), m_ionic_param(size) {} //! \name Size //! \brief Manage the size of the list // -------- //! @{ //! \brief Resize the list void resize(index_t size) override { ReactiveSpeciesList::resize(size); m_ionic_param.resize(size); } //! \brief Resize the list, adjust the number of components void resize(index_t size, index_t nb_components) override { ReactiveSpeciesList::resize(size, nb_components); m_ionic_param.resize(size); } //! @} //! \name Getter //! \brief Return values // ---------------------- //! @{ //! \brief Return the ionic size of component 'k' const scalar_t& a_debye(index_t k) const { return m_ionic_param.a_debye(k); } //! \brief Return the 'b-dot' parameter of component 'k' const scalar_t& b_debye(index_t k) const { return m_ionic_param.b_debye(k); } //! \brief Return the charge of component 'k' const scalar_t& charge(index_t k) const { return m_ionic_param.charge(k); } //! \brief Return the ionic model values of component 'k' IonicModelValues ionic_values(index_t k) const { return m_ionic_param.get_values(k); } //! @} //! \name Setter //! \brief Set values // ------------- //! @{ //! \brief Set the ionic model parameters void set_ionic_values(index_t k, const IonicModelValues& values) { m_ionic_param.set_values(k, values); } //! \brief Set the values //! \warning Do no set stoichiometric coefficients void set_values(index_t k, const AqueousValues& values) { set_label(k, values.label); set_logk(k, values.logk); m_ionic_param.set_values(k, values.ionic_values); } //! \brief Set the values //! \warning Do no set stoichiometric coefficients void set_values(index_t k, AqueousValues&& values) { set_label(k, std::move(values.label)); set_logk(k, values.logk); m_ionic_param.set_values(k, std::move(values.ionic_values)); } //! @} //! \name Move //! \brief Move species inside the list or to other lists // ------------ //! @{ //! \brief Move component 'old_ind' to 'new_ind' void move_erase(index_t old_ind, index_t new_ind) override { m_ionic_param.move_erase(old_ind, new_ind); ReactiveSpeciesList::move_erase(old_ind, new_ind); } //! \brief Move component 'old_ind' to 'new_ind' void move_erase(index_t old_ind, index_t new_ind, const std::vector& is_reactants_to_remove) override { m_ionic_param.move_erase(old_ind, new_ind); ReactiveSpeciesList::move_erase(old_ind, new_ind, is_reactants_to_remove); } //! @} //! \brief Add the aqueous species 'other_species' to 'k', to obtain a canonical system void canonicalize(index_t k, index_t other_species, scalar_t coeff) { add_species_to(k, other_species, coeff); } + + friend class specmicp::AdimensionalSystem; + private: IonicModelParameters m_ionic_param; }; } // end namespace database } // end namespace specmicp #endif //SPECMICP_DATABASE_AQUEOUSLIST_HPP diff --git a/src/specmicp_database/species/base_wrapper.hpp b/src/specmicp_database/species/base_wrapper.hpp index 06d6dba..e12f859 100644 --- a/src/specmicp_database/species/base_wrapper.hpp +++ b/src/specmicp_database/species/base_wrapper.hpp @@ -1,269 +1,275 @@ /* ============================================================================= 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_DATABASE_BASEWRAPPER_HPP #define SPECMICP_DATABASE_BASEWRAPPER_HPP #include "specmicp_common/types.hpp" #include //! \file base_wrapper.hpp //! //! \brief Base wrapper around vectors and matrices //! //! These classes encapsulate an Eigen Matrix (or vector) //! //! \defgroup database_species namespace specmicp { + +class AdimensionalSystem; + namespace database { //! \class MatrixSpeciesWrapper //! \brief A wrapper around a matrix, to be used by a SpeciesList (or a derived class) instance //! //! \ingroup database_species Describing the species in the database class SPECMICP_DLL_PUBLIC MatrixSpeciesWrapper { public: //! \brief Default constructor MatrixSpeciesWrapper() {} //! \brief Constructor initializing the matrix MatrixSpeciesWrapper(index_t size, index_t nb_cols): m_matrix(Matrix::Zero(size, nb_cols)) {} //! \name Size // ============ //! \brief Return and set the size of the matrix //! @{ //! \brief Return the number of rows in the matrix index_t size() const {return m_matrix.rows();} //! \brief Return the number of columnns in the matrix index_t cols() const {return m_matrix.cols();} //! \brief Resize the matrix //! This is a resize only for the number of rows void resize(index_t size) { m_matrix.conservativeResize(size, Eigen::NoChange); } //! \brief Resize the matrix (both rows and colums) void resize(index_t rows, index_t cols) { m_matrix.conservativeResize(rows, cols); } // @} //! \name Getter // ============= //! \brief Return the values //! @{ //! \brief Return the element ['row', ['col'] const scalar_t& operator()(index_t row, index_t col) const { return m_matrix(row, col); } //! \brief Return a const reference to the matrix const Matrix& get_matrix() const { return m_matrix; } //! \brief Return the row of a matrix Eigen::MatrixBase::ConstRowXpr get_row(index_t k) const { return m_matrix.row(k); } //! \brief Return the element ['row','col'] const scalar_t& get_value(index_t row, index_t col) const { return m_matrix(row, col); } //! @} //! \name Setter // ============= //! \brief Set the values //! @{ //! \brief Add 'other' to row 'k' template void SPECMICP_DLL_LOCAL add_to_row(index_t k, const Eigen::MatrixBase& other) { m_matrix.row(k) += other; } //! \brief Multiply row 'k' by 'multiplier' template void SPECMICP_DLL_LOCAL multiply_by(const Eigen::MatrixBase& multiplier) { m_matrix = m_matrix*multiplier; } //! \brief Set the value of the matrix by copying 'init_matrix' template void SPECMICP_DLL_LOCAL set_matrix(const Eigen::MatrixBase& init_matrix) { m_matrix = init_matrix; } //! \brief Set the element ['row','col'] to 'value' void SPECMICP_DLL_LOCAL set_value(index_t row, index_t col, scalar_t value) { m_matrix(row, col) = value; } //! \brief Reset row 'k' void SPECMICP_DLL_LOCAL reset_row(index_t k) { m_matrix.row(k).setZero(); } //! @} //! \name Move // =========== //! \brief Move values inside, and outside, of the matrix //! @{ //! \brief replace imformation of 'new_ind' by those contained in 'old_ind' void move_erase(index_t old_ind, index_t new_ind) { m_matrix.row(new_ind) = m_matrix.row(old_ind); } //! \brief replace imformation of 'new_ind' by those contained in 'old_ind' void move_erase(index_t old_ind, index_t new_ind, const std::vector& is_col_to_remove); //! \brief Move the row to another matrix void move_erase_to(index_t ind, MatrixSpeciesWrapper& other, index_t other_ind) { other.m_matrix.row(other_ind) = m_matrix.row(ind); m_matrix.row(ind).setZero(); } //! \brief Swap two rows void swap(index_t ind, MatrixSpeciesWrapper& other, index_t other_ind) { m_matrix.row(ind).swap(other.m_matrix.row(other_ind)); } //! @} + //! + + friend class specmicp::AdimensionalSystem; private: Matrix m_matrix; }; //! \class VectorSpeciesWrapper //! \brief A wrapper around a vector //! //! \ingroup database_species class SPECMICP_DLL_PUBLIC VectorSpeciesWrapper { public: //! \brief Default constructor VectorSpeciesWrapper() {} //! \brief Initialise the vector //! \param size initial size of the vector VectorSpeciesWrapper(index_t size): m_vector(Vector::Zero(size)) {} //! \name Size // ============ //! \brief Return and set the size of the vector //! @{ //! \brief Return the size of the vector index_t size() const {return m_vector.rows();} //! \brief Resize the vector void resize(index_t size) { m_vector.conservativeResize(size); } //! @} //! \name Getter // ============= //! \brief Return the values //! @{ //! \brief return the value of index k const scalar_t& operator()(index_t k) const { return m_vector(k); } //! \brief Return the inner product of the vector and 'other' template scalar_t dot(const Eigen::MatrixBase& other) const { return m_vector.dot(other); } //! \brief Return the value of index k const scalar_t& get_value(index_t k) const { return m_vector(k); } //! \brief Return the vector const Vector& get_vector() const { return m_vector; } //! @} //! \name Setter // ============= //! \brief Set the values //! @{ //! \brief add 'vector' template void SPECMICP_DLL_LOCAL add(const Eigen::MatrixBase& vector) { m_vector += vector; } //! \brief Multiply the vector by 'multiplier' template void SPECMICP_DLL_LOCAL multiply_by(const Eigen::MatrixBase& multiplier) { m_vector = m_vector*multiplier; } //! \brief Set the value at index k void SPECMICP_DLL_LOCAL set_value(index_t k, scalar_t value) { m_vector(k) = value; } //! \brief Copy 'vector' to the vector template void SPECMICP_DLL_LOCAL set_vector(const Eigen::MatrixBase& vector) { m_vector = vector; } //! \brief Apply a linear transformation to the vector template void transform(const Eigen::MatrixBase& transform_matrix) { specmicp_assert(transform_matrix.rows() == transform_matrix.cols()); m_vector = transform_matrix*m_vector; } //! @} //! \name Move // =========== //! \brief Move values inside, and outside, of the matrix //! @{ //! \brief Move value of old_ind at new_ind void move_erase(index_t old_ind, index_t new_ind) { m_vector(new_ind) = m_vector(old_ind); } //! \brief Move row 'ind' to 'other' at row 'other_ind' void move_erase_to(index_t ind, VectorSpeciesWrapper& other, index_t other_ind) { other.m_vector(other_ind) = m_vector(ind); m_vector(ind) = 0; } //! @} private: Vector m_vector; }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_BASEWRAPPER_HPP diff --git a/src/specmicp_database/species/component_list.hpp b/src/specmicp_database/species/component_list.hpp index 53c033d..e1dbde5 100644 --- a/src/specmicp_database/species/component_list.hpp +++ b/src/specmicp_database/species/component_list.hpp @@ -1,165 +1,170 @@ /* ============================================================================= 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_DATABASE_COMPONENTLIST_HPP #define SPECMICP_DATABASE_COMPONENTLIST_HPP #include "specmicp_common/types.hpp" #ifndef SPECMICP_DATABASE_SPECIES_HPP #include "species.hpp" #endif #ifndef SPECMICP_DATABASE_IONICPARAMETERS_HPP #include "ionic_parameters.hpp" #endif //! \file component_list.hpp A list of components, i.e the basis namespace specmicp { + +class AdimensionalSystem; + namespace database { //! \struct ComponentValues //! \brief Used for the initialization of a component //! //! \ingroup database_species struct ComponentValues { std::string label; scalar_t molar_mass; IonicModelValues ionic_values; }; //! \class ComponentList //! \brief The basis //! //! This is the list of components in the system //! It extends the species list by managing the ionic //! model parameters and the molar masses. //! //! \ingroup database_species class ComponentList: public SpeciesList { public: ComponentList() {} ComponentList(index_t size): SpeciesList(size), m_molar_mass(size), m_ionic_param(size) {} // Resize // ------ //! \brief Resize the basis void resize(index_t size) override { SpeciesList::resize(size); m_molar_mass.resize(size); m_ionic_param.resize(size); } // Getter // ====== //! \brief Return the molar mass of component k const scalar_t& molar_mass(index_t k) const { return m_molar_mass(k); } //! \brief Return the charge of component 'k' const scalar_t& charge(index_t k) const { return m_ionic_param.charge(k); } //! \brief Return the ionic size of component 'k' const scalar_t& a_debye(index_t k) const { return m_ionic_param.a_debye(k); } //! \brief Return the 'b-dot' parameter of component 'k' const scalar_t& b_debye(index_t k) const { return m_ionic_param.b_debye(k); } //! \brief Return the ionic model values of component 'k' IonicModelValues ionic_values(index_t k) const { return m_ionic_param.get_values(k); } // Setter // ------ //! \brief set the values for component 'k' void set_values(index_t k, const ComponentValues& values) { set_label(k, values.label); m_molar_mass.set_value(k, values.molar_mass); m_ionic_param.set_values(k, values.ionic_values); } //! \brief set the values for component 'k' void set_values(index_t k, ComponentValues&& values) { set_label(k, std::move(values.label)); m_molar_mass.set_value(k, values.molar_mass); m_ionic_param.set_values(k, values.ionic_values); } void set_ionic_values(index_t k, const IonicModelValues& values) { m_ionic_param.set_values(k, values); } //! \brief Return the molar mass of a compounts //! \param stoech a vector containing the stoichiometric coefficients template scalar_t molar_mass_compounds(const Eigen::MatrixBase& stoech) const { return m_molar_mass.dot(stoech); } //! \brief Transform the molar mass of the basis //! \param transform_matrix the linear operator to apply template void molar_mass_transform(const Eigen::MatrixBase& transform_matrix) { m_molar_mass.transform(transform_matrix); } // Move // ----- //! \brief Move component 'old_ind' to 'new_ind' void move_erase(index_t old_ind, index_t new_ind) override { m_molar_mass.move_erase(old_ind, new_ind); m_ionic_param.move_erase(old_ind, new_ind); SpeciesList::move_erase(old_ind, new_ind); } + friend class specmicp::AdimensionalSystem; + private: VectorSpeciesWrapper m_molar_mass; IonicModelParameters m_ionic_param; }; } // end namespace database } // end namespace specmicp #endif //SPECMICP_DATABASE_COMPONENTLIST_HPP diff --git a/src/specmicp_database/species/ionic_parameters.hpp b/src/specmicp_database/species/ionic_parameters.hpp index a85b73d..10a9dcf 100644 --- a/src/specmicp_database/species/ionic_parameters.hpp +++ b/src/specmicp_database/species/ionic_parameters.hpp @@ -1,114 +1,120 @@ /* ============================================================================= 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_DATABASE_IONICPARAMETERS_HPP #define SPECMICP_DATABASE_IONICPARAMETERS_HPP #ifndef SPECMICP_DATABASE_BASEWRAPPER_HPP #include "base_wrapper.hpp" #endif //! \file ionic_parameters.hpp //! \brief The ionic model parameters //! //! These parameters are needed for all (primary and secondary) //! aqueous species. namespace specmicp { + +class AdimensionalSystem; + namespace database { //! \struct IonicModelValues //! The set of main values for the ionic model //! This is used for the initialisation of an aqueous species. //! //! \ingroup database_species struct IonicModelValues { scalar_t charge {0.0}; scalar_t a_debye {0.0}; scalar_t b_debye {0.0}; IonicModelValues() {} IonicModelValues(scalar_t acharge, scalar_t adebye, scalar_t bdebye): charge(acharge), a_debye(adebye), b_debye(bdebye) {} }; //! \class IonicModelParameters //! \brief Class that contains the ionic model parameters //! //! This class is used for all (primary and secondary) aqueous species. //! //! \ingroup database_species class IonicModelParameters: public MatrixSpeciesWrapper { public: static constexpr index_t charge_ind = 0; //!< Column index of the charge static constexpr index_t adebye_ind = 1; //!< Column index of the ionic size parameter static constexpr index_t bdebye_ind = 2; //!< Column index of the b-dot parameter IonicModelParameters() {} IonicModelParameters(index_t size): MatrixSpeciesWrapper(size, 3) {} void resize_columns(index_t colums) = delete; //!< The number of columns is fixed void resize_matrix(index_t rows, index_t cols) = delete; //!< The number of columns is fixed //! \brief Return the charge const scalar_t& charge(index_t ind) const {return get_value(ind, charge_ind);} //! \brief Return the ionic size parameter const scalar_t& a_debye(index_t ind) const {return get_value(ind, adebye_ind);} //! \brief Return the b-dot parameter const scalar_t& b_debye(index_t ind) const {return get_value(ind, bdebye_ind);} //! \brief Set the values void set_values(index_t ind, const IonicModelValues& values) { set_value(ind, charge_ind, values.charge); set_value(ind, adebye_ind, values.a_debye); set_value(ind, bdebye_ind, values.b_debye); } //! \brief Return the ionic value IonicModelValues get_values(index_t ind) const { return IonicModelValues(charge(ind), a_debye(ind), b_debye(ind)); } + + + friend class specmicp::AdimensionalSystem; }; } // end namespace database } // end namespace specmicp #endif //SPECMICP_DATABASE_IONICPARAMETERS_HPP