diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index ffd99ed..e103f6b 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,1553 +1,1553 @@ /* ============================================================================= 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; } } else if (constraints.water_equation == WaterEquationType::FixedSaturation) { m_fixed_values(dof_water()) = constraints.water_parameter; } else { specmicp_assert(constraints.water_equation == WaterEquationType::SaturatedSystem); } m_equations.add_equation(dof_water(), &neq); } // Aqueous components // ================== number_eq_aqueous_component(constraints, neq); // Surface model // ============= 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; // true if all the components of the mineral // are present in the system // just check that the molar volume exist auto molar_volume = m_data->molar_volume_mineral(m); // Remove minerals that cannot precipitate for (index_t& k: m_equations.nonactive_component) { if (m_data->nu_mineral(m, k) != 0.0 and molar_volume > 0.0) { can_precipitate = false; break; // this is not a mineral that can precipitate } } // check if no precipitation is set // if upper bound is 0, deactivate the mineral if (not constraints.mineral_constraints.empty()) { for (const MineralConstraint& v: constraints.mineral_constraints) { if (v.id_mineral == m) { if (v.param == 0) { can_precipitate = false; } } } } if (can_precipitate) { m_equations.add_equation(dof_mineral(m), &neq); } } m_equations.nb_tot_variables = neq; m_equations.nb_complementarity_variables = m_equations.nb_tot_variables - m_equations.nb_free_variables; // Minerals => no precipitation // // /!\ must be set after equations ids are set if (not constraints.mineral_constraints.empty()) { for (const MineralConstraint& v: constraints.mineral_constraints) { if (v.equation_type == MineralConstraintType::NoPrecipitation) { index_t idm = v.id_mineral; if (idm == no_species or ideq_min(idm) == no_equation) continue; if (not m_equations.is_box_constrained) { m_equations.activate_box_constrained(neq); } m_equations.set_box_constrained(ideq_min(idm), v.param); } } } } void AdimensionalSystem::number_eq_aqueous_component( const AdimensionalSystemConstraints& constraints, index_t& neq ) { using EqT = AqueousComponentEquationType; // First set the charge keeper if (constraints.charge_keeper != no_species) { if (constraints.charge_keeper == 0 or constraints.charge_keeper > m_data->nb_component()) { throw std::invalid_argument("The charge keeper must be an aqueous component. Invalid argument : " + std::to_string(constraints.charge_keeper)); } m_equations.type_equation(dof_component(constraints.charge_keeper)) = static_cast(EqT::ChargeBalance); } // Then go over fix fugacity gas for (const auto& it: constraints.fixed_fugacity_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component) + "' is already constrained, a fixed fugacity condition can not be applied"); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedFugacity); m_fixed_values(it.id_component) = it.log_value; m_equations.related_species(it.id_component) = it.id_gas; } // Then over the fixed activity species for (const auto& it: constraints.fixed_activity_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component) + "' is already constrained, a fixed activity condition can not be applied."); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedActivity); m_fixed_values(it.id_component) = it.log_value; } // Then the fixed molality components for (const auto& it: constraints.fixed_molality_cs) { if (m_equations.type_equation(dof_component(it.id_component)) != static_cast(EqT::NoEquation)) { throw std::invalid_argument("Component '" + m_data->components.get_label(it.id_component) + "' is already constrained, a fixed molality condition can not be applied."); } m_equations.type_equation(dof_component(it.id_component)) = static_cast(EqT::FixedMolality); m_fixed_values(it.id_component) = it.log_value; } // 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) + // (latter is necessary because units are always reset by solver) m_scaling_molar_volume = m_data->scaling_molar_volume(get_units()); // Unit scaling for the gaseous total concentration // transform mol/m^3 into the right unit (pressure are always in Pa) switch (get_units().length) { case units::LengthUnit::decimeter: m_scaling_gas = 1e-3; break; case units::LengthUnit::centimeter: m_scaling_gas = 1e-6; break; default: m_scaling_gas = 1.0; break; } if (get_units().quantity == units::QuantityUnit::millimoles) { m_scaling_gas *= 1e3; m_scaling_molality = 1e3; } } // ================ // // // // 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 aporosity = m_second.porosity; scalar_t sat_w = volume_fraction_water(x) / aporosity; if (sat_w < 1) { const scalar_t pressure = m_equations.water_pressure_model(sat_w); res -= m_scaling_gas*(aporosity - volume_fraction_water(x))*( pressure / (constants::gas_constant*temperature())); } } res /= total_concentration_bc(0); return res; } scalar_t AdimensionalSystem::residual_water_saturation(const Vector& x) const { specmicp_assert(water_equation_type() == WaterEquationType::SaturatedSystem); scalar_t res = 1 - volume_fraction_water(x) - m_inert_volume_fraction; for (index_t mineral: m_data->range_mineral()) { res -= volume_fraction_mineral(x, mineral); } return res; } scalar_t AdimensionalSystem::residual_water_fixed_saturation(const Vector& x) const { specmicp_assert(water_equation_type() == WaterEquationType::FixedSaturation); scalar_t porosity = 1 - m_inert_volume_fraction; for (index_t mineral: m_data->range_mineral()) { porosity -= volume_fraction_mineral(x, mineral); } scalar_t res = fixed_saturation_bc()*porosity - volume_fraction_water(x); return res; } scalar_t AdimensionalSystem::residual_component(const Vector &x, index_t component) const { specmicp_assert(aqueous_component_equation_type(component) == AqueousComponentEquationType::MassConservation); const scalar_t conc_w = density_water()*volume_fraction_water(x); scalar_t res = total_concentration_bc(component); scalar_t aqconc = component_molality(x, component) + weigthed_sum_aqueous(component); if (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 = porosity(x); set_volume_fraction_gas_phase(x); set_secondary_concentration(x); set_sorbed_concentrations(x); set_pressure_fugacity(x); // // water if (ideq_w() != no_equation) { switch (water_equation_type()) { case WaterEquationType::MassConservation: residual(ideq_w()) = residual_water_conservation(x); break; case WaterEquationType::SaturatedSystem: residual(ideq_w()) = residual_water_saturation(x); break; case WaterEquationType::FixedSaturation: residual(ideq_w()) = residual_water_fixed_saturation(x); break; case WaterEquationType::NoEquation: break; } } // aqueous component for (index_t i: m_data->range_aqueous_component()) { switch (aqueous_component_equation_type(i)) { case AqueousComponentEquationType::NoEquation: break; case AqueousComponentEquationType::MassConservation: residual(ideq_paq(i)) = residual_component(x, i); break; case AqueousComponentEquationType::ChargeBalance: residual(ideq_paq(i)) = residual_charge_conservation(x); break; case AqueousComponentEquationType::FixedActivity: residual(ideq_paq(i)) = residual_fixed_activity(x, i); break; case AqueousComponentEquationType::FixedFugacity: residual(ideq_paq(i)) = residual_fixed_fugacity(x, i); break; case AqueousComponentEquationType::FixedMolality: residual(ideq_paq(i)) = residual_fixed_molality(x, i); break; } } // surface if (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; } } else if (water_equation_type() == WaterEquationType::FixedSaturation) { const index_t idw = ideq_w(); jacobian(idw, idw) = -1; for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idw, ideq_min(m)) = - fixed_saturation_bc(); } } } void AdimensionalSystem::jacobian_aqueous_components(Vector& x, Matrix& jacobian) { for (index_t i: m_data->range_aqueous_component()) { const index_t idp = ideq_paq(i); if (idp == no_equation) continue; switch (aqueous_component_equation_type(i)) { case AqueousComponentEquationType::NoEquation: continue; // Mass balance equation // ===================== case AqueousComponentEquationType::MassConservation: { const scalar_t conc_w = density_water()*volume_fraction_water(x); const scalar_t factor = total_concentration_bc(i); // Aqueous components for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp_iip = 0; if (k == i) tmp_iip -= component_molality(x, i)*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()); vector_debye_huckel_component(ion_str, sqrti); vector_debye_huckel_aqueous(ion_str, sqrti); } void AdimensionalSystem::vector_debye_huckel_component( scalar_t ionic_strength, scalar_t sqrt_ionic ) { Vector& log_gamma = m_second.loggamma; Matrix& ionic_param = m_data->components.m_ionic_param.m_matrix; using IoParam = database::IonicModelParameters; for (auto i: m_data->range_aqueous_component()) { const scalar_t zisquare = std::pow(ionic_param(i, IoParam::charge_ind), 2); log_gamma(i) = - ( constants::Adebye * zisquare * sqrt_ionic); log_gamma(i) /= ( 1 + ionic_param(i, IoParam::adebye_ind)*constants::Bdebye*sqrt_ionic); log_gamma(i) += ionic_param(i, IoParam::bdebye_ind)*ionic_strength; } } void AdimensionalSystem::vector_debye_huckel_aqueous( scalar_t ionic_strength, scalar_t sqrt_ionic ) { Vector& log_gamma = m_second.loggamma; Matrix& ionic_param = m_data->aqueous.m_ionic_param.m_matrix; using IoParam = database::IonicModelParameters; const index_t offset = m_data->nb_component(); for (auto j: m_data->range_aqueous()) { const scalar_t zisquare = std::pow(ionic_param(j, IoParam::charge_ind), 2); log_gamma(offset+j) = - (constants::Adebye * zisquare * sqrt_ionic); log_gamma(offset+j) /= ( 1 + ionic_param(j, IoParam::adebye_ind)*constants::Bdebye*sqrt_ionic); log_gamma(offset+j) += ionic_param(j, IoParam::bdebye_ind)*ionic_strength; } } bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual) { if (not get_options().non_ideality) { set_secondary_variables(x); // we still need to compute secondary species ! // if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) // { // set_secondary_variables(x); // } return true; } not_in_linesearch = true; scalar_t previous_norm = m_second.loggamma.norm(); if (previous_norm == 0) previous_norm = 1.0; bool may_have_converged = false; if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) { m_equations.solve_pressure_model = true; // Use fixed point iterations for non-ideality for (int i=0; i 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) { static std::mt19937 gen(std::random_device{}()); std::uniform_real_distribution<> dis(-2, 2); xtot(dof_water()) = get_options().restart_water_volume_fraction; for (index_t i: m_data->range_aqueous_component()) { //if (xtot(dof_component(i)) > 0 or xtot(dof_component(i)) < -9) xtot(i) = get_options().restart_concentration + dis(gen); } if (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_common/log.hpp b/src/specmicp_common/log.hpp index 0c3f082..4dd114f 100644 --- a/src/specmicp_common/log.hpp +++ b/src/specmicp_common/log.hpp @@ -1,164 +1,164 @@ /* ============================================================================= 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_UTILS_LOG_HPP #define SPECMICP_UTILS_LOG_HPP #include /*! \file log.hpp \brief logger A common logger is shared by all SpecMiCP/ReactMiCP module. The logger contains several level which can be activated (or not) by the user. The user also decides of where the logger output its information. \code{.cpp} // output to std::cerr for Warning and Errors (nmo Debug messages) init_logger(&std::cerr, logger::LogLevel::Warning) // Log an Error ERROR << "Error"; // Log a Warning WARNING << "Warning"; // Log a debug message (no effect); DEBUG << "Debug (disabled message) \endcode This module was inspired by : - http://www.drdobbs.com/cpp/a-lightweight-logger-for-c/240147505 - http://www.drdobbs.com/cpp/logging-in-c/201804215?pgno=1 */ #include "log_impl.hpp" namespace specmicp { //! \namespace specmicp::logger //! \brief The logger used by SpecMiCP/ReactMiCP namespace logger { extern template class Log; extern template class Log; } // end namespace logger //! \brief Standard logger type using stdlog = logger::Log; using conflog = logger::Log; //! \brief Initialize the logger //! //! \param out the output stream //! \param level the output level void init_logger(std::ostream* out, specmicp::logger::LogLevel level); //! \brief log a message and throw a runtime error void log_and_throw_runtime_error(const std::string& error_msg); - //! \brief Configure the special logger containing information about //! the configuration of specmicp void init_conf_logger(std::ostream* out); + #define SPC_CONF_LOG_FILTER(level) \ if (logger::ConfFile::stream() == nullptr) ;\ else conflog().get(level) #define SPC_CONF_LOG SPC_CONF_LOG_FILTER(logger::Info) #define SPC_CONF_LOG_WARNING SPC_CONF_LOG_FILTER(logger::Warning) #define SPC_CONF_LOG_ERROR SPC_CONF_LOG_FILTER(logger::Error) #define SPC_CONF_LOG_CRITICAL SPC_CONF_LOG_FILTER(logger::Critical) #define SPC_CONF_LOG_SECTION "======================================" #define SPC_CONF_LOG_HLINE "--------------------------------------" //! \brief Filter logs to stdlog #define SPC_LOG_FILTER(level) \ if (level > stdlog::ReportLevel() || logger::ErrFile::stream() == nullptr) ;\ else stdlog().get(level) #define SPC_LOG_FILTER_NDEBUG \ if (true) ; \ else stdlog().get(logger::Warning) //! \def SPAM //! \brief Output spam-level log to stdlog #ifdef NDEBUG #define SPAM SPC_LOG_FILTER_NDEBUG #else #define SPAM SPC_LOG_FILTER(logger::Spam) #endif //! \def DEBUG //! \brief Output debug-level log to stdlog #ifdef NDEBUG #define DEBUG SPC_LOG_FILTER_NDEBUG #else #define DEBUG SPC_LOG_FILTER(logger::Debug) #endif //! \def INFO //! \brief Ooutput info-level log to stdlog #ifdef INFO #undef INFO // jsoncpp is also using this, remove compiler warnings #endif #ifdef NDEBUG #define INFO SPC_LOG_FILTER_NDEBUG #else #undef INFO #define INFO SPC_LOG_FILTER(logger::Info) #endif //! \def WARNING //! \brief Output warning-level log to stdlog #define WARNING SPC_LOG_FILTER(logger::Warning) //! \def ERROR //! \brief Output error-level log to stdlog #define ERROR SPC_LOG_FILTER(logger::Error) //! \def ERROR_THROW //! \brief Output error level log to stdlog and throw runtime error #define ERROR_THROW(msg) log_and_throw_runtime_error(msg) //! \def CRITICAL //! \brief Output critical-level log to stdlog #define CRITICAL SPC_LOG_FILTER(logger::Critical) } // end namespace specmicp #endif // SPECMICP_UTILS_LOG_HPP diff --git a/src/specmicp_common/physics/unsaturated_laws.hpp b/src/specmicp_common/physics/unsaturated_laws.hpp index 97e3e76..2e50c55 100644 --- a/src/specmicp_common/physics/unsaturated_laws.hpp +++ b/src/specmicp_common/physics/unsaturated_laws.hpp @@ -1,132 +1,132 @@ /* ============================================================================= 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_PHYSICS_UNSATURATEDLAWS_HPP #define SPECMICP_PHYSICS_UNSATURATEDLAWS_HPP //! \file physics/unsaturated_laws.hpp //! \brief Common laws and models in unsaturated porous medium #include "specmicp_common/types.hpp" namespace specmicp { namespace laws { //! \brief Return the capillary pressure as function of the relative humidity scalar_t SPECMICP_DLL_PUBLIC kelvin_equation(scalar_t rh); //! \brief Return the relative humidity as function of the capillary pressure scalar_t SPECMICP_DLL_PUBLIC invert_kelvin_equation(scalar_t pc); //! \brief Capillary Pressure model used by V. Baroghel-Bouny et al. //! //! \f$ Pc(S_l) = a (S_l^{-b} - 1)^{1-1/b} \f$ //! //! \param saturation the liquid saturation //! \param a fitting coefficient, unit of pressure //! \param b fitting coefficient, dimensioneless, >1 scalar_t SPECMICP_DLL_PUBLIC capillary_pressure_BaroghelBouny( scalar_t saturation, scalar_t a, scalar_t b ); //! \brief Van Genutchen capillary Pressure model //! //! \f$ Pc(S_l) = P_0 (S_l^{-1/m} - 1)^{1-m} \f$ //! //! \param saturation the liquid saturation //! \param P_0 fitting coefficient, unit of pressure //! \param m fitting coefficient, dimensioneless, 0; -// FIXME don't work with GCC < 4.9.0 +// don't work with GCC < 4.9.0 / no fix //! \brief Return a model function of the saturation //! //! \tparam F a model where saturation is the first parameter //! \tparam ...Args Parameter pack representing the parameters of the model //! \return a function taking only the saturation //! //! Bind the argument of the model template saturation_model_f make_saturation_model(F func, Args...args) { return [func,args...](scalar_t saturation) -> scalar_t { return func(saturation, args...); }; } //! \brief Return a relative humidity model from a capillary pressure model //! //! Use the kelvin equation inline saturation_model_f make_rh_model(saturation_model_f cap_pressure) { return [cap_pressure](scalar_t saturation) -> scalar_t { return invert_kelvin_equation(cap_pressure(saturation)); }; } } //end namespace laws } //end namespace specmicp #endif // SPECMICP_PHYSICS_UNSATURATEDLAWS_HPP diff --git a/src/specmicp_database/yaml_reader.cpp b/src/specmicp_database/yaml_reader.cpp index 3ebfd99..27bdba1 100644 --- a/src/specmicp_database/yaml_reader.cpp +++ b/src/specmicp_database/yaml_reader.cpp @@ -1,749 +1,750 @@ /* ============================================================================= 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 "config_database.hpp" #include "yaml_reader.hpp" #include "reader_common.hpp" #include "section_name.hpp" #include "specmicp_common/io/yaml.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/compat.hpp" #include #include #include namespace specmicp { namespace database { struct DataReaderYaml::DataReaderYamlElementData { bool check_compo {true}; std::vector components; bool is_init {false}; }; DataReaderYaml::DataReaderYaml(): m_elem_data(make_unique()) {} DataReaderYaml::~DataReaderYaml() = default; //! \brief Constructor DataReaderYaml::DataReaderYaml(RawDatabasePtr data, bool check_compo): DatabaseModule(data), m_elem_data(make_unique()) { m_elem_data->check_compo = check_compo; } //! \brief Constructor //! //! @param filepath string containing the path to the database DataReaderYaml::DataReaderYaml(const std::string& filepath, bool check_compo): DatabaseModule(), m_elem_data(make_unique()) { m_elem_data->check_compo = check_compo; parse(filepath); } //! \brief Constructor //! //! @param input input stream that contains the database DataReaderYaml::DataReaderYaml(std::istream& input, bool check_compo): DatabaseModule(), m_elem_data(make_unique()) { m_elem_data->check_compo = check_compo; parse(input); } // safe access to values std::string obtain_label(const YAML::Node& species, const std::string& section); void obtain_composition(const YAML::Node& species, std::map& compo, const std::string& label); scalar_t obtain_logk(const YAML::Node& species, const std::string& label); bool is_kinetic(const YAML::Node& species,const std::string& label); void DataReaderYaml::parse(const std::string& filepath) { std::ifstream datafile(filepath, std::ifstream::binary); if (not datafile.is_open()) { ERROR << "Database not found : " << filepath; std::string message("Database not found : " + filepath); throw std::invalid_argument(message); } data->metadata.path = filepath; parse(datafile); datafile.close(); } void DataReaderYaml::parse(std::istream& input ) { const YAML::Node root = io::parse_yaml(input); const char* indb_main = INDB_MAIN; io::check_mandatory_yaml_node(root, INDB_SECTION_METADATA, indb_main); parse_metadata(root[INDB_SECTION_METADATA]); // Basis // ------ io::check_mandatory_yaml_node(root, INDB_SECTION_BASIS, indb_main); parse_basis(root[INDB_SECTION_BASIS]); // Elements // -------- io::check_mandatory_yaml_node(root, INDB_SECTION_ELEMENTS, indb_main); ElementList elements; parse_elements(root[INDB_SECTION_ELEMENTS], elements); data->elements = std::move(elements); // Aqueous species // --------------- io::check_mandatory_yaml_node(root, INDB_SECTION_AQUEOUS, indb_main); AqueousList alist; parse_aqueous(root[INDB_SECTION_AQUEOUS], alist); data->aqueous = std::move(alist); // Minerals // -------- io::check_mandatory_yaml_node(root, INDB_SECTION_MINERALS, indb_main); MineralList minerals, minerals_kinetic; parse_minerals(root[INDB_SECTION_MINERALS], minerals, minerals_kinetic); data->minerals = std::move(minerals); data->minerals_kinetic = std::move(minerals_kinetic); // Gas // ---- GasList glist(0, data->nb_component()); if (root[INDB_SECTION_GAS]) parse_gas(root[INDB_SECTION_GAS], glist); else glist.set_valid(); data->gas = std::move(glist); // Sorbed species // -------------- SorbedList slist(0, data->nb_component()); if (root[INDB_SECTION_SORBED]) parse_sorbed(root[INDB_SECTION_SORBED], slist); else slist.set_valid(); data->sorbed = std::move(slist); // Compounds // --------- CompoundList clist(0, data->nb_component()); if (root[INDB_SECTION_COMPOUNDS]) parse_compounds(root[INDB_SECTION_COMPOUNDS], clist); else clist.set_valid(); data->compounds = std::move(clist); } void DataReaderYaml::parse_metadata(const YAML::Node& root) { - data->metadata.name = - io::get_yaml_mandatory(root, INDB_ATTRIBUTE_NAME, INDB_SECTION_METADATA); - data->metadata.version = - io::get_yaml_mandatory(root, INDB_ATTRIBUTE_VERSION, INDB_SECTION_METADATA); - - m_elem_data->check_compo = io::get_yaml_optional(root, INDB_ATTRIBUTE_CHECKCOMPO, INDB_SECTION_METADATA, m_elem_data->check_compo); + data->metadata.name = io::get_yaml_mandatory( + root, INDB_ATTRIBUTE_NAME, INDB_SECTION_METADATA); + data->metadata.version = io::get_yaml_mandatory( + root, INDB_ATTRIBUTE_VERSION, INDB_SECTION_METADATA); + m_elem_data->check_compo = io::get_yaml_optional( + root, INDB_ATTRIBUTE_CHECKCOMPO, INDB_SECTION_METADATA, + m_elem_data->check_compo); } void DataReaderYaml::init_elemental_composition() { if (m_elem_data->is_init) return; m_elem_data->components.reserve(data->nb_component()); // Parse label to find elemental composition for (auto id: data->range_component()) { element_map elem_compo; element_composition_from_label(data->get_label_component(id), elem_compo); m_elem_data->components.emplace_back(elem_compo); } m_elem_data->is_init = true; } void DataReaderYaml::parse_basis(const YAML::Node& basis) { DEBUG << "Size basis : " << basis.size(); data->components = ComponentList(basis.size()); if (data->nb_component() <= 3) throw std::invalid_argument("The basis is too small."); index_t starting_point = 2; for (int id: data->components.range()) { index_t id_in_db = starting_point; if (id_in_db >= data->nb_component()) { throw std::invalid_argument("The basis should contain H2O and E[-]."); } const YAML::Node& species = basis[id]; std::string label = obtain_label(species, INDB_SECTION_BASIS); if (label == water_label) { id_in_db = 0; //throw invalid_database("The first component should be "+water_label+"."); } else if (label == electron_label) { id_in_db = 1; //throw invalid_database("The second component should be "+electron_label+"."); } else { ++starting_point; } auto is_already_in_the_database = data->components.get_id(label); if (is_already_in_the_database != no_species) throw invalid_database("Component : " + label + "is already in the database."); SPAM << "Parsing component : " << label; ComponentValues values {}; values.label = label; // activity coeeficients values.ionic_values.charge = charge_from_label(label); if (species[INDB_ATTRIBUTE_ACTIVITY]) { values.ionic_values.a_debye = io::get_yaml_mandatory( species[INDB_ATTRIBUTE_ACTIVITY], INDB_ATTRIBUTE_ACTIVITY_A, label); values.ionic_values.b_debye = io::get_yaml_mandatory( species[INDB_ATTRIBUTE_ACTIVITY], INDB_ATTRIBUTE_ACTIVITY_B, label); } else { values.ionic_values.a_debye = 0; values.ionic_values.b_debye = 0; } // molar mass values.molar_mass = io::get_yaml_mandatory( species, INDB_ATTRIBUTE_MOLARMASS, label); data->components.set_values(id_in_db, std::move(values)); } if (m_elem_data->check_compo) init_elemental_composition(); // after this the basis is ready to use data->components.set_valid(); } void DataReaderYaml::parse_aqueous(const YAML::Node& aqueous, AqueousList& alist) { DEBUG << "Parse aqueous species"; alist = AqueousList(aqueous.size(), data->nb_component()); //const auto size = data->nb_aqueous(); for (auto id: alist.range()) { const YAML::Node& species = aqueous[static_cast(id)]; AqueousValues values; values.label = obtain_label(species, INDB_SECTION_AQUEOUS); values.ionic_values.charge = charge_from_label(values.label); if (species[INDB_ATTRIBUTE_ACTIVITY]) { values.ionic_values.a_debye = io::get_yaml_mandatory( species[INDB_ATTRIBUTE_ACTIVITY], INDB_ATTRIBUTE_ACTIVITY_A, values.label); values.ionic_values.b_debye = io::get_yaml_mandatory( species[INDB_ATTRIBUTE_ACTIVITY], INDB_ATTRIBUTE_ACTIVITY_B, values.label); } values.logk = obtain_logk(species, values.label); alist.set_values(id, values); // equation std::map compo; obtain_composition(species, compo, values.label); scalar_t test_charge {0.0}; std::map to_canonicalize; for (const auto& it: compo) { // component ? const auto search = data->components.get_id(it.first); if(search != no_species) { alist.set_nu_ji(id, search, it.second); test_charge += it.second*data->charge_component(search); } // aqueous species ? else { // in the document we parse ? const auto id_aq = alist.get_id(it.first); if (id_aq != no_species) { to_canonicalize.insert({id_aq, it.second}); test_charge += it.second*alist.charge(id_aq); } // cannot be found... else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + alist.get_label(id) +"."); } } } if (std::abs(test_charge - alist.charge(id)) > EPS_TEST_CHARGE) { throw db_invalid_data("Total charge is not zero for aqueous species : '" + alist.get_label(id) + "' - Charge-decomposition - charge species : "+ std::to_string(test_charge - alist.charge(id) )+"."); } for (const auto& tocanon: to_canonicalize) { alist.canonicalize(id, tocanon.first, tocanon.second); } // check composition if (m_elem_data->check_compo) { if (species[INDB_ATTRIBUTE_FORMULA]) check_composition(species[INDB_ATTRIBUTE_FORMULA].as(), alist, id); else check_composition(alist.get_label(id), alist, id); } } // valid set of aqueous species alist.set_valid(); } void DataReaderYaml::parse_minerals( const YAML::Node& minerals, MineralList& minerals_list, MineralList& minerals_kinetic_list ) { // this one is a little bit more complicated since we need to detect // if it is a solid phase governed by equilibrium or kinetic DEBUG << "Parse minerals"; index_t nb_mineral = minerals.size(); // find kinetic minerals int nb_kin = 0; for (int id=0; idnb_component()); minerals_kinetic_list = MineralList(nb_kin, data->nb_component()); // id for each class of minerals index_t id_in_eq=0; index_t id_in_kin=0; for (index_t id=0; id(id)]; //check if material is at equilibrium or governed by kinetics MineralValue value; value.label = obtain_label(species, INDB_SECTION_MINERALS); value.logk = obtain_logk(species, value.label); bool is_kin = is_kinetic(species, value.label); // ###FIXME value.molar_volume = io::get_yaml_optional(species, INDB_ATTRIBUTE_MOLARVOLUME, value.label, -1); auto& mlist = is_kin?minerals_kinetic_list:minerals_list; auto& true_id = is_kin?id_in_kin:id_in_eq; mlist.set_values(true_id, value); // equation // -------- std::map compo; obtain_composition(species, compo, value.label); double test_charge = 0; std::map to_canonicalize; for (const auto& it: compo) { auto idc = data->components.get_id(it.first); if(idc != no_species) { // this is a component mlist.set_nu_ji(true_id, idc, it.second); test_charge += it.second * data->charge_component(idc); } else { // this is an aqueous species auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second * data->charge_aqueous(idaq); } else // or something we don't know { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + mlist.get_label(true_id) +"."); } } } if (std::abs(test_charge) > EPS_TEST_CHARGE) { throw db_invalid_data("Total charge is not zero for mineral : '"+ mlist.get_label(true_id) + "' - total charge : "+std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { mlist.canonicalize(true_id, data->aqueous, tocanon.first, tocanon.second); } // check composition if (m_elem_data->check_compo) { if (species[INDB_ATTRIBUTE_FORMULA]) check_composition(species[INDB_ATTRIBUTE_FORMULA].as(), mlist, true_id); else check_composition(mlist.get_label(true_id), mlist, true_id); } // update index ++true_id; } specmicp_assert(id_in_kin == minerals_kinetic_list.size()); specmicp_assert(id_in_eq == minerals_list.size()); // ok after that point minerals_list.set_valid(); minerals_kinetic_list.set_valid(); } void DataReaderYaml::parse_gas(const YAML::Node& gas, GasList& glist) { DEBUG << "Parse gas"; glist = GasList(gas.size(), data->nb_component()); for (auto id: glist.range()) { const YAML::Node& species = gas[static_cast(id)]; GasValues values; values.label = obtain_label(species, INDB_SECTION_GAS); auto is_already_in_the_database = data->gas.get_id(values.label); if (is_already_in_the_database != no_species) throw invalid_database("Component : " + values.label + "is already in the database."); values.logk = obtain_logk(species, values.label); glist.set_values(id, values); // equation std::map compo; obtain_composition(species, compo, values.label); scalar_t test_charge = 0.0; std::map to_canonicalize; for (auto it: compo) { const auto idc = data->components.get_id(it.first); // component if(idc != no_species) { glist.set_nu_ji(id, idc, it.second); test_charge += it.second*data->charge_component(idc); } // aqueous species else { const auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second*data->charge_aqueous(idaq); } else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + data->gas.get_label(id) +"."); } } } if (std::abs(test_charge) > EPS_TEST_CHARGE) { throw db_invalid_data("Total charge is not zero for gas '"+data->gas.get_label(id)+ "' : " + std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { glist.canonicalize(id, data->aqueous, tocanon.first, tocanon.second); } // check composition if (m_elem_data->check_compo) { if (species[INDB_ATTRIBUTE_FORMULA]) check_composition(species[INDB_ATTRIBUTE_FORMULA].as(), glist, id); else check_composition(glist.get_label(id), glist, id); } } glist.set_valid(); } void DataReaderYaml::parse_sorbed(const YAML::Node& sorbed, SorbedList& slist) { DEBUG << "Parse sorbed species"; slist = SorbedList(sorbed.size(), data->nb_component()); for (auto id: slist.range()) { const YAML::Node& species = sorbed[static_cast(id)]; SorbedValues values; values.label = obtain_label(species, INDB_SECTION_SORBED); auto is_already_in_the_database = slist.get_id(values.label); if (is_already_in_the_database != no_species) throw invalid_database("Sorbed species : " + values.label + "is already in the database."); values.logk = obtain_logk(species, values.label); const scalar_t nb_site_occupied = io::get_yaml_mandatory (species, INDB_ATTRIBUTE_NBSITEOCCUPIED, values.label); if (nb_site_occupied < 0) { throw db_invalid_data("The number of sites occupied by a sorbed species must be positive. (Species : " + values.label + ", number of sites : " +std::to_string(nb_site_occupied) + ")"); } values.sorption_site_occupied = nb_site_occupied; slist.set_values(id, values); std::map compo; obtain_composition(species, compo, values.label); double test_charge = 0; std::map to_canonicalize; for (auto it: compo) { const auto idc = data->components.get_id(it.first); if(idc != no_species) { slist.set_nu_ji(id, idc, it.second); test_charge += it.second*data->charge_component(idc); } else { const auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second*data->charge_aqueous(idaq); } else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + slist.get_label(id) +"."); } } } if (std::abs(test_charge) > EPS_TEST_CHARGE) { throw db_invalid_data("Total charge is not zero for gas '"+slist.get_label(id)+ "' : " + std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { slist.canonicalize(id, data->aqueous, tocanon.first, tocanon.second); } } slist.set_valid(); } void DataReaderYaml::parse_compounds(const YAML::Node& compounds, CompoundList &clist) { DEBUG << "Parse compounds"; clist = CompoundList(compounds.size(), data->nb_component()); for (auto id: clist.range()) { const YAML::Node& species = compounds[static_cast(id)]; std::string label = obtain_label(species, INDB_SECTION_COMPOUNDS); auto is_already_in_the_database = clist.get_id(label); if (is_already_in_the_database != no_species) throw invalid_database("Compounds : " + label + "is already in the database."); clist.set_values(id, label); std::map compo; obtain_composition(species, compo, label); double test_charge = 0; std::map to_canonicalize; for (auto it: compo) { const auto idc = data->components.get_id(it.first); if(idc != no_species) { clist.set_nu_ji(id, idc, it.second); test_charge += it.second*data->charge_component(idc); } else { const auto idaq = data->aqueous.get_id(it.first); if (idaq != no_species) { to_canonicalize.insert({idaq, it.second}); test_charge += it.second*data->charge_aqueous(idaq); } else { throw db_invalid_syntax("Unknown species : '" + it.first + "' while parsing equations for " + clist.get_label(id) +"."); } } } if (std::abs(test_charge) > EPS_TEST_CHARGE) { throw db_invalid_data("Total charge is not zero for gas '"+clist.get_label(id)+ "' : " + std::to_string(test_charge)+"."); } // canonicalise for (const auto& tocanon: to_canonicalize) { clist.canonicalize(id, data->aqueous, tocanon.first, tocanon.second); } // check composition if (m_elem_data->check_compo) { if (species[INDB_ATTRIBUTE_FORMULA]) check_composition(species[INDB_ATTRIBUTE_FORMULA].as(), clist, id); else check_composition(clist.get_label(id), clist, id); } } clist.set_valid(); } void check_element_map(const element_map& formula, const element_map& composition, const std::string& label); void DataReaderYaml::check_composition(const std::string& formula, ReactiveSpeciesList& rlist, index_t id) { if (not m_elem_data->is_init) init_elemental_composition(); // parse the formula element_map compo_from_label; element_composition_from_label(formula, compo_from_label); // parse the composition element_map compo_from_compo; for (auto idc: data->range_component()) { if (rlist.nu_ji(id, idc) == 0.0) continue; add_to_element_map(compo_from_compo, m_elem_data->components[idc], rlist.nu_ji(id, idc) ); } // clean the composition for (auto it=compo_from_compo.begin(); it!=compo_from_compo.end(); ) { if (it->first == "E" or it->second == 0.0) it = compo_from_compo.erase(it); else ++it; } // check that both map are equal check_element_map(compo_from_label, compo_from_compo, rlist.get_label(id)); } void check_element_map(const element_map& formula, const element_map& composition, const std::string& label) { for (const auto& it1: formula) { const auto it2 = composition.find(it1.first); if (it2 == composition.cend()) { throw db_invalid_data("Species '"+label+"' : element '"+ it1.first+"' is in the formula but not in the composition"); } if (std::abs(it1.second - it2->second) > EPS_TEST_COMPOSITION) { throw db_invalid_data("Species '"+label+"' : element '"+ it1.first+"' has a different stoichiometry in the formula ("+ std::to_string(it1.second)+ ") than in the composition ("+ std::to_string(it2->second)+")."); } } for (const auto& it1: composition) { const auto it2 = formula.find(it1.first); if (it2 == composition.cend()) { throw db_invalid_data("Species '"+label+"' : element '"+ it1.first+"' is in the composition but not in the formula"); } } } void DataReaderYaml::parse_elements(const YAML::Node& elements, ElementList &elist) { if (elements.size() != static_cast(data->nb_component())) throw db_invalid_data("The number of elements does not corresponds to the number of components"); for (index_t id: data->range_component()) { const YAML::Node& elem = elements[static_cast(id)]; const std::string elem_label = io::get_yaml_mandatory( elem, INDB_ATTRIBUTE_ELEMENT, INDB_SECTION_ELEMENTS); const std::string comp_label = io::get_yaml_mandatory( elem, INDB_ATTRIBUTE_COMPONENT, elem_label); index_t id_comp = data->get_id_component(comp_label); if (id_comp == no_species) { throw db_invalid_label("'"+comp_label+"' is not a valid component"); } elist.add_element(elem_label, id_comp); } } std::string obtain_label(const YAML::Node& species, const std::string& section) { return io::get_yaml_mandatory(species, INDB_ATTRIBUTE_LABEL, section); } void obtain_composition(const YAML::Node& species, std::map& compo, const std::string& label) { const std::string compo_string = io::get_yaml_mandatory(species, INDB_ATTRIBUTE_COMPOSITION, label); parse_equation(compo_string, compo); } scalar_t obtain_logk(const YAML::Node& species, const std::string& label) { return io::get_yaml_mandatory(species, INDB_ATTRIBUTE_LOGK, label); } bool is_kinetic(const YAML::Node& species, const std::string& label) { return io::get_yaml_optional(species, INDB_ATTRIBUTE_FLAG_KINETIC, label, false); } } //end namespace database } //end namespace specmicp