diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index 874718f..6913b22 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,1430 +1,1430 @@ /* ============================================================================= 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() { m_scaling_molar_volume = m_data->scaling_molar_volume(get_units().length); // 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; } } // ================ // // // // 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); res = total_concentration_bc(0); res -= conc_w/m_data->molar_mass_basis(0); res -= conc_w*weigthed_sum_aqueous(0); if (ideq_surf() != no_equation) res -= conc_w*weigthed_sum_sorbed(0); res -= weigthed_sum_mineral(x, 0); if (m_data->nb_gas() > 0) res -= weigthed_sum_gas(0); if (m_equations.use_water_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); res -= conc_w*component_molality(x, component); res -= conc_w*weigthed_sum_aqueous(component); if (ideq_surf() != no_equation) res -= conc_w*weigthed_sum_sorbed(component); 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 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 -= 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); // // 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 *= rho_w; tmp -= -weigthed_sum_gas(0); if (m_equations.use_water_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) { // doesn't make sense to have a saturation // bigger than one in this case ERROR << "Saturation greater that one detected" << ", skip pressure model"; } 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 -= diff_weigthed_sum_aqueous(k, 0); tmp -= 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 -= diff_weigthed_sum_aqueous(m_data->electron_index(), 0); tmp -= 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 *= 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 *= density_water(); jacobian(idp, ideq_w()) = tmp_iw/factor; } // Surface if (ideq_surf() != no_equation) { scalar_t tmp_s = -conc_w*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*= 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) += 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) = 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) = 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; set_volume_fraction_gas_phase(x); set_pressure_fugacity(x); if (ideq_surf() != no_equation) set_sorbed_concentrations(x); set_secondary_concentration(x); compute_log_gamma(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 sqrti = std::sqrt(ionic_strength()); for (index_t 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) ); } for (index_t 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) ); } } 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; bool may_have_converged = false; if (norm_residual < nb_free_variables()*get_options().start_non_ideality_computation) { // 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(std::abs(previous_norm - m_second.loggamma.norm())); } } // Set the correct value for the water total saturation if (ideq_w() == no_equation) { xtot(dof_water()) = volume_fraction_water(x); } return AdimensionalSystemSolution(xtot, m_second.secondary_molalities, m_second.loggamma, m_second.ionic_strength, m_second.gas_fugacity, m_second.sorbed_concentrations, m_inert_volume_fraction); } // Water, saturation and density // ============================== scalar_t AdimensionalSystem::density_water() const { return laws::density_water(units::celsius(25.0), length_unit(), mass_unit()); } 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()) = 0.8; + xtot(dof_water()) = get_options().restart_water_volume_fraction; for (index_t i: m_data->range_aqueous_component()) { xtot(dof_component(i)) = -4.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()) = 0.5; + 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_structs.hpp b/src/specmicp/adimensional/adimensional_system_structs.hpp index 130c6dc..161ce1d 100644 --- a/src/specmicp/adimensional/adimensional_system_structs.hpp +++ b/src/specmicp/adimensional/adimensional_system_structs.hpp @@ -1,340 +1,342 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP #include "specmicp_common/types.hpp" //! \file adimensional_system_structs.hpp //! \brief Options and constraints for the AdimensionalSystem namespace specmicp { //! \struct AdimensionalSystemOptions //! \brief Options for the Adimensional Systems //! //! It is mainly about the secondary variables fixed-point iterations struct SPECMICP_DLL_PUBLIC AdimensionalSystemOptions { //! Solve for non ideality bool non_ideality; //! Max iterations for the non ideality model index_t non_ideality_max_iter; //! Scaling the electron equation scalar_t scaling_electron; //! Tolerance for non ideality scalar_t non_ideality_tolerance; //! Under relaxation factor for the conservation of water scalar_t under_relaxation_factor; //! Log of the molality used to restart the computation scalar_t restart_concentration; + //! Liquid water volume fraction to restart the computation + scalar_t restart_water_volume_fraction; //! Log_10 of the molality for a new component scalar_t new_component_concentration; //! Factor to start the non-ideality computation scalar_t start_non_ideality_computation; //! Cutoff for including components in the computation scalar_t cutoff_total_concentration; AdimensionalSystemOptions(); // in adimensional_system_solver_structs.cpp }; //! \enum AqueousComponentEquationType //! \brief Type of an aqueous component equation enum class AqueousComponentEquationType { NoEquation = no_equation, //!< Not an equation, component is not present in the system MassConservation, //!< Mass balance ChargeBalance, //!< M.B. replaced by charge balance FixedFugacity, //!< M.B. replaced by a fixed fugacity equation FixedActivity, //!< M.B. replaced by a fixed activity equation FixedMolality //!< The molality of the species is fixed }; //! \enum WaterEquationType //! \brief The type of the equation solved for the water enum class WaterEquationType { NoEquation = no_equation, //!< Amount of water is not solved MassConservation, //!< Water is conserved SaturatedSystem //!< System is saturated }; //! \struct FixedFugacityConstraint //! \brief Struct to contain information needed to solve a fix fugacity problem struct FixedFugacityConstraint { index_t id_gas; //!< Index of the fixed-fugacity gas index_t id_component; //!< Index of the corresponding component scalar_t log_value; //!< Log_10 of the fugacity FixedFugacityConstraint(index_t gas, index_t component, scalar_t logvalue): id_gas(gas), id_component(component), log_value(logvalue) {} }; //! \struct FixedActivityConstraint //! \brief Struct to contain information needed to solve a fix activity problem. struct FixedActivityConstraint { index_t id_component; //!< Index of the fixed-activity component scalar_t log_value; //!< Log_10 of the activity FixedActivityConstraint(index_t component, scalar_t logvalue): id_component(component), log_value(logvalue) {} }; //! \struct FixedMolalityConstraint //! \brief Contain information about a fixed molality constraint struct FixedMolalityConstraint { index_t id_component; //!< Index of the component scalar_t log_value; //!< Log_10 of the molality FixedMolalityConstraint(index_t component, scalar_t logvalue): id_component(component), log_value(logvalue) {} }; //! \enum SurfaceEquationType //! \brief The model for surface sorption enum class SurfaceEquationType { NoEquation = no_equation, //!< Do not include surface sorption Equilibrium //!< Equilibrium model }; //! \struct SurfaceConstraint //! This struct contains the information to set-up the surface sorption model struct SurfaceConstraint { SurfaceEquationType model_type; //!< The model to use scalar_t concentration; //!< The total concentration of sorption sites //! \brief By default, we don't include surface sorption in the computation SurfaceConstraint() noexcept: model_type(SurfaceEquationType::NoEquation), concentration(0.0) {} //! \brief When a concentration is supplied, the surface sorption model is equilibrium SurfaceConstraint(scalar_t surface_concentration) noexcept: model_type(SurfaceEquationType::Equilibrium), concentration(surface_concentration) {} }; //! \enum ElectronEquationType the type of the equation for the electron enum class ElectronEquationType { NoEquation = no_equation, //!< Do not compute the concentration equation of the electron Equilibrium, //!< Set the concentration of electron to be 0 FixedpE //!< Activity of the electron is fixed }; //! \struct ElectronConstraint //! \brief the constraint for the electron struct ElectronConstraint { ElectronEquationType equation_type{ElectronEquationType::NoEquation}; //!< The equation type scalar_t fixed_value; //!< The fixed value of pE if needed index_t species {no_species}; //!< In case of fixed pE, this is the reaction to use //! \brief By default we assume equilibrium ElectronConstraint(): equation_type(ElectronEquationType::Equilibrium), fixed_value(0.0) {} //! \brief When a value is provided, we assume that the pE is fixed ElectronConstraint(scalar_t pe_value, scalar_t aqueous_species): equation_type(ElectronEquationType::FixedpE), fixed_value(pe_value), species(aqueous_species) {} }; //! \brief Return the partial pressure as function of the liquid saturation //! //! This function must return the pressure in Pa using water_partial_pressure_f = std::function; //! \brief Constraints for the partial pressure model struct WaterPartialPressureConstraint { bool use_partial_pressure_model; //!< True if we use partial pressure model water_partial_pressure_f partial_pressure_model; //!< The model given by the user //! \brief By default the partial pressure constraint is disabled WaterPartialPressureConstraint(): use_partial_pressure_model(false), partial_pressure_model(nullptr) {} //! \brief Enble the water partial pressure model WaterPartialPressureConstraint(water_partial_pressure_f& model): WaterPartialPressureConstraint() { if (model != nullptr) { use_partial_pressure_model = true; partial_pressure_model = model; } } }; //! \struct AdimensionalSystemConstraints //! \brief Struct to contains the "Boundary conditions" for the AdimensionalSystem //! //! \ingroup specmicp_api struct AdimensionalSystemConstraints { Vector total_concentrations; //!< Total concentrations WaterEquationType water_equation{WaterEquationType::MassConservation}; //!< Water equation index_t charge_keeper{no_species}; //!< The equation for this component is replace by the charge balance std::vector fixed_fugacity_cs {}; //!< Contains information about fixed fugacity constraints std::vector fixed_activity_cs {}; //!< Contains information about fixed activity constraints std::vector fixed_molality_cs {}; //!< Contains information about fixed molality constraints scalar_t inert_volume_fraction{ 0.0 }; //!< Volume fraction of inert solid (inert in the equilibrium computation) SurfaceConstraint surface_model{}; //!< Surface sorption model ElectronConstraint electron_constraint{}; //!< constraint for the electron WaterPartialPressureConstraint water_partial_pressure {}; //!< The partial pressure of water (if any) AdimensionalSystemConstraints(): total_concentrations() {} AdimensionalSystemConstraints(const Vector& total_concs): total_concentrations(total_concs) {} //! \brief Set the total concentrations void set_total_concentrations(const Vector& total_concs) {total_concentrations = total_concs;} //! \brief Enable the conservation of water void enable_conservation_water() noexcept {water_equation = WaterEquationType::MassConservation;} //! \brief Disable the conservation of water void disable_conservation_water() noexcept {water_equation = WaterEquationType::NoEquation;} //! \brief The system is saturated void set_saturated_system() noexcept {water_equation = WaterEquationType::SaturatedSystem;} //! \brief Disable the surface sorption model void disable_surface_model() noexcept {surface_model.model_type = SurfaceEquationType::NoEquation;} //! \brief Enable the surface sorption model //! \param surface_sorption_model_concentration concentration of the surface sorption sites void enable_surface_model(scalar_t surface_sorption_model_concentration) noexcept { surface_model.model_type = SurfaceEquationType::Equilibrium; surface_model.concentration = surface_sorption_model_concentration; } //! \brief Set the charge keeper to 'component' //! //! \param component Index of the component (in the database) void set_charge_keeper(index_t component) noexcept { charge_keeper = component; } //! \brief Add a fixed fugacity gas condition //! //! \param constraint struct containing the information about a fixed-fugacity constraint void add_fixed_fugacity_gas(const FixedFugacityConstraint& constraint) { fixed_fugacity_cs.push_back(constraint); } //! \brief Add a fixed fugacity gas condition //! //! \param gas Index of the gas (in the database) //! \param component Index of the corresponding component (in the database) //! \param logvalue Log_10 of the fugacity void add_fixed_fugacity_gas(index_t gas, index_t component, scalar_t logvalue) { fixed_fugacity_cs.emplace_back(FixedFugacityConstraint(gas, component, logvalue)); } //! \brief Add a fixed activity component condition //! //! \param constraint struct containing the information about a fixed-activity constraint void add_fixed_activity_component(const FixedActivityConstraint& constraint) { fixed_activity_cs.push_back(constraint); } //! \brief Add a fixed activity condition for a component //! //! \param component Index of the corresponding component (in the database) //! \param log_value Log_10 of the activity void add_fixed_activity_component(index_t component, scalar_t log_value) { fixed_activity_cs.emplace_back(FixedActivityConstraint(component, log_value)); } //! \brief Add a fixed molality condition for a component //! //! \param component Index of the corresponding component (in the database) //! \param log_value Log_10 of the molality void add_fixed_molality_component(index_t component, scalar_t log_value) { fixed_molality_cs.emplace_back(FixedMolalityConstraint(component, log_value)); } //! \brief Set the inert volume fraction //! //! The volume fraction of the inert phase is used to offset the saturation. //! This inert phase may correspond to aggregates or solid phases governed by kinetics. //! //! \param value volume fraction of the inert phase void set_inert_volume_fraction(scalar_t value) noexcept { inert_volume_fraction = value; } //! \brief Set a partial pressure model for water //! //! The partial pressure model for water is a function taking the saturation //! and returning the partial pressure of water void set_water_partial_pressure_model(water_partial_pressure_f model) { water_partial_pressure.use_partial_pressure_model = true; water_partial_pressure.partial_pressure_model = model; } // //! \brief Set the system at a fixed pE // void set_fixed_pe(scalar_t pe_value, index_t aqueous_species) noexcept { // electron_constraint = ElectronConstraint(pe_value, aqueous_species); // } }; } // end namespace specmicp #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP diff --git a/src/specmicp/adimensional/config_default_options_solver.h b/src/specmicp/adimensional/config_default_options_solver.h index 3205591..b0ae9b8 100644 --- a/src/specmicp/adimensional/config_default_options_solver.h +++ b/src/specmicp/adimensional/config_default_options_solver.h @@ -1,89 +1,90 @@ /* ============================================================================= 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. * ============================================================================= */ //! \file config_default_options_solver.h //! \internal //! \brief Defaut options for the adimensional solver #ifndef SPECMICP_ADIM_CONFIG_DEFAULTOPTIONSSOLVER_H #define SPECMICP_ADIM_CONFIG_DEFAULTOPTIONSSOLVER_H #include "specmicp_common/micpsolver/micpsolver_structs.hpp" // System options #define SPECMICP_DEFAULT_ENABLE_NONIDEAL true #define SPECMICP_DEFAULT_NONIDEAL_MAX_ITER 10 #define SPECMICP_DEFAULT_SCALING_ELECTRON 0.0 #define SPECMICP_DEFAULT_NONIDEAL_TOL 1e-8 #define SPECMICP_DEFAULT_UNDERRELAX_FACTOR 0.9 #define SPECMICP_DEFAULT_RESTART_CONC -6.0 +#define SPECMICP_DEFAULT_RESTART_WATER 0.5 #define SPECMICP_DEFAULT_NEW_COMPONENT_CONC -4.0 #define SPECMICP_DEFAULT_TRHSOLD_START_NONIDEAL 0.1 #define SPECMICP_DEFAULT_CUTOFF_TOT_CONC 1e-12 // Micpslover options #define SPECMICP_DEFAULT_USE_CRASHING MICPSOLVER_DEFAULT_USE_CRASHING #define SPECMICP_DEFAULT_USE_SCALING MICPSOLVER_DEFAULT_USE_SCALING #define SPECMICP_DEFAULT_USE_NONMONOTONE_LSEARCH MICPSOLVER_DEFAULT_USE_NONMONOTONE_LSEARCH #define SPECMICP_DEFAULT_MAX_ITER MICPSOLVER_DEFAULT_MAX_ITER #define SPECMICP_DEFAULT_MAX_ITER_MAXSTEP MICPSOLVER_DEFAULT_MAX_ITER_MAXSTEP #define SPECMICP_DEFAULT_MAX_FACT_STEP MICPSOLVER_DEFAULT_MAX_FACT_STEP #define SPECMICP_DEFAULT_RES_TOL MICPSOLVER_DEFAULT_RES_TOL #define SPECMICP_DEFAULT_STEP_TOL MICPSOLVER_DEFAULT_STEP_TOL #define SPECMICP_DEFAULT_COND_THRESHOLD MICPSOLVER_DEFAULT_COND_THRESHOLD #define SPECMICP_DEFAULT_PENALIZATION_FACTOR MICPSOLVER_DEFAULT_PENALIZATION_FACTOR #define SPECMICP_DEFAULT_MAX_STEP 10.0 #define SPECMICP_DEFAULT_FACTOR_DESC_COND -1.0 #define SPECMICP_DEFAULT_POWER_DESC_COND MICPSOLVER_DEFAULT_POWER_DESC_COND #define SPECMICP_DEFAULT_COEFF_ACCEPT_NEWTON MICPSOLVER_DEFAULT_COEFF_ACCEPT_NEWTON #define SPECMICP_DEFAULT_FACTOR_GRADIENT MICPSOLVER_DEFAULT_FACTOR_GRADIENT #define SPECMICP_DEFAULT_PROJ_VAR MICPSOLVER_DEFAULT_PROJ_VAR #define SPECMICP_DEFAULT_TRHSOLD_STATIONARY MICPSOLVER_DEFAULT_TRHSOLD_STATIONARY #define SPECMICP_DEFAULT_TRHSOLD_CYCLING_LSEARCH MICPSOLVER_DEFAULT_TRHSOLD_CYCLING_LSEARCH // Adim solver options #define SPECMICP_DEFAULT_ALLOW_RESTART true #define SPECMICP_DEFAULT_USE_PCFM false #define SPECMICP_DEFAULT_FORCE_PCFM false #define SPECMICP_DEFAULT_LENGTH_UNIT specmicp::units::LengthUnit::meter #endif // SPECMICP_ADIM_CONFIG_DEFAULTOPTIONSSOLVER_H diff --git a/src/specmicp/io/configuration.cpp b/src/specmicp/io/configuration.cpp index 370ab56..9b666e7 100644 --- a/src/specmicp/io/configuration.cpp +++ b/src/specmicp/io/configuration.cpp @@ -1,898 +1,900 @@ /* ============================================================================= 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 "configuration.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/adimensional/config_default_options_solver.h" #include "specmicp_common/io/yaml.hpp" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp_common/physics/laws.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/io/config_yaml_sections.h" namespace specmicp { namespace io { // Additional declarations // ======================= void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_activity, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_fugacity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_fugacity, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_molality( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_molality, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_activity ); void configure_specmicp_constraints_fixed_fugacity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_fugacity ); void configure_specmicp_constraints_fixed_molality( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_molality ); scalar_t molar_mass_aqueous_unknown_class(const std::string& label, const RawDatabasePtr& raw_db); scalar_t molar_mass_solid_unknown_class(const std::string& label, const RawDatabasePtr& raw_db ); scalar_t molar_volume_solid_unknown_class(const std::string& label, const RawDatabasePtr& raw_db ); // Implementation // ============== void configure_specmicp_options( AdimensionalSystemSolverOptions& options, const YAML::Node& configuration ) { check_mandatory_yaml_node(configuration, SPC_CF_S_SPECMICP, "__main__"); const YAML::Node& conf = configuration[SPC_CF_S_SPECMICP]; // fixme units //options.units_set = the_units_set; options.solver_options.fvectol = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_RES_TOL, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_RES_TOL); options.solver_options.steptol = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_STEP_TOL, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_STEP_TOL); options.solver_options.max_iter = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_MAX_ITER, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_MAX_ITER); options.solver_options.maxstep = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_MAX_STEP); options.solver_options.maxiter_maxstep = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER, SPC_CF_S_SPECMICP, std::min(SPECMICP_DEFAULT_MAX_ITER_MAXSTEP, options.solver_options.max_iter) ); options.solver_options.use_scaling = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_SCALING, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_USE_SCALING); options.solver_options.non_monotone_linesearch = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_NONMONOTONE, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_USE_NONMONOTONE_LSEARCH); options.solver_options.factor_descent_condition = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_FACTOR_DESC_COND ); options.solver_options.condition_limit = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_COND_CHECK, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_COND_THRESHOLD); options.solver_options.threshold_cycling_linesearch = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_TRHSOLD_CYCLING_LSEARCH ); options.system_options.non_ideality_tolerance = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_NONIDEAL_TOL, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_NONIDEAL_TOL); options.system_options.cutoff_total_concentration = get_yaml_optional(conf, SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC, SPC_CF_S_SPECMICP, SPECMICP_DEFAULT_CUTOFF_TOT_CONC); if (configuration[SPC_CF_S_UNITS]) { options.units_set = configure_units(configuration[SPC_CF_S_UNITS]); } } void configure_specmicp_options( AdimensionalSystemSolverOptions& options, const units::UnitsSet& the_units, YAMLConfigHandle&& conf ) { micpsolver::MiCPSolverOptions& solver = options.solver_options; conf.set_if_attribute_exists( solver.fvectol, SPC_CF_S_SPECMICP_A_RES_TOL, 0.0); conf.set_if_attribute_exists( solver.steptol, SPC_CF_S_SPECMICP_A_STEP_TOL, 0.0); conf.set_if_attribute_exists( solver.max_iter, SPC_CF_S_SPECMICP_A_MAX_ITER, 0); conf.set_if_attribute_exists( solver.maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH, 0.0); conf.set_if_attribute_exists( solver.maxiter_maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER, 0); conf.set_if_attribute_exists( solver.use_scaling, SPC_CF_S_SPECMICP_A_SCALING); conf.set_if_attribute_exists( solver.non_monotone_linesearch, SPC_CF_S_SPECMICP_A_NONMONOTONE); conf.set_if_attribute_exists( solver.factor_descent_condition, SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION); conf.set_if_attribute_exists( solver.condition_limit, SPC_CF_S_SPECMICP_A_COND_CHECK); conf.set_if_attribute_exists( solver.threshold_cycling_linesearch, SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH); AdimensionalSystemOptions& system = options.system_options; conf.set_if_attribute_exists( system.non_ideality_tolerance, SPC_CF_S_SPECMICP_A_NONIDEAL_TOL, 0.0); conf.set_if_attribute_exists( system.non_ideality_max_iter, SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER, 0); conf.set_if_attribute_exists( system.cutoff_total_concentration, SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC, 0.0); conf.set_if_attribute_exists( system.restart_concentration, SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION); + conf.set_if_attribute_exists( + system.restart_water_volume_fraction, SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC, 0.0); conf.set_if_attribute_exists( system.under_relaxation_factor, SPC_CF_S_SPECMICP_A_UNDER_RELAXATION, 0.0, 1.0); options.units_set = the_units; } void configure_specmicp_formulation( Formulation& formulation, const YAML::Node& conf_formulation, const RawDatabasePtr& raw_db, const units::UnitsSet& units_sets ) { configure_specmicp_formulation_solution(formulation, conf_formulation, units_sets); configure_specmicp_formulation_aqueous(formulation, conf_formulation, raw_db, units_sets); configure_specmicp_formulation_solid(formulation, conf_formulation, raw_db, units_sets); } void configure_specmicp_formulation_solution( Formulation& formulation, const YAML::Node& conf_formulation, const units::UnitsSet& _) { check_mandatory_yaml_node(conf_formulation, SPC_CF_S_FORMULATION_A_SOLUTION, SPC_CF_S_FORMULATION); const YAML::Node& sol_conf = conf_formulation[SPC_CF_S_FORMULATION_A_SOLUTION]; check_mandatory_yaml_node(sol_conf, SPC_CF_S_FORMULATION_A_AMOUNT, SPC_CF_S_FORMULATION_A_SOLUTION); scalar_t factor=1.0; if (sol_conf[SPC_CF_S_FORMULATION_A_UNIT]) { units::AmountUnit unit; units::parse_amount_unit(sol_conf[SPC_CF_S_FORMULATION_A_UNIT].as(), unit); if (unit.type != units::AmountUnitType::Mass ) { throw std::invalid_argument("'" + sol_conf[SPC_CF_S_FORMULATION_A_UNIT].as() + "'' is not a mass unit (in the solution configuration) !"); } factor = unit.factor_si; } formulation.set_mass_solution(sol_conf[SPC_CF_S_FORMULATION_A_AMOUNT].as()*factor); } void configure_specmicp_formulation_aqueous( Formulation& formulation, const YAML::Node& conf_formulation, const RawDatabasePtr& raw_db, const units::UnitsSet& _) { check_mandatory_yaml_node( conf_formulation, SPC_CF_S_FORMULATION_A_AQUEOUS, SPC_CF_S_FORMULATION); const YAML::Node& aq_conf = conf_formulation[SPC_CF_S_FORMULATION_A_AQUEOUS]; for (auto aqueous_conf: aq_conf) { check_mandatory_yaml_node(aqueous_conf, SPC_CF_S_FORMULATION_A_LABEL, SPC_CF_S_FORMULATION_A_AQUEOUS); check_mandatory_yaml_node(aqueous_conf, SPC_CF_S_FORMULATION_A_AMOUNT, SPC_CF_S_FORMULATION_A_AQUEOUS); check_mandatory_yaml_node(aqueous_conf, SPC_CF_S_FORMULATION_A_UNIT, SPC_CF_S_FORMULATION_A_AQUEOUS); std::string label = aqueous_conf[SPC_CF_S_FORMULATION_A_LABEL].as(); scalar_t factor = 1.0; units::AmountUnit unit; units::parse_amount_unit(aqueous_conf[SPC_CF_S_FORMULATION_A_UNIT].as(), unit); if (unit.type == units::AmountUnitType::Molality) { factor = unit.factor_si; } else if (unit.type == units::AmountUnitType::MoleConcentration) { factor = unit.factor_si/laws::density_water(units::celsius(25.0), units::LengthUnit::meter, units::MassUnit::kilogram); } else if (unit.type == units::AmountUnitType::MassConcentration) { const scalar_t molar_mass = molar_mass_aqueous_unknown_class(label, raw_db); factor = unit.factor_si / molar_mass / laws::density_water(units::celsius(25.0), units::LengthUnit::meter, units::MassUnit::kilogram); } else if (unit.type == units::AmountUnitType::Mass) { const scalar_t molar_mass = molar_mass_aqueous_unknown_class(label, raw_db); factor = unit.factor_si / molar_mass / formulation.mass_solution; } else if (unit.type == units::AmountUnitType::NumberOfMoles) { factor = unit.factor_si / formulation.mass_solution; } else { throw std::invalid_argument("'"+aqueous_conf[SPC_CF_S_FORMULATION_A_UNIT].as() + "' is not a valid unit for an aqueous species"); } const scalar_t molality = factor*aqueous_conf[SPC_CF_S_FORMULATION_A_AMOUNT].as(); formulation.add_aqueous_species(label, molality); } } void configure_specmicp_formulation_solid( Formulation& formulation, const YAML::Node& conf_formulation, const RawDatabasePtr& raw_db, const units::UnitsSet& unit_set) { check_mandatory_yaml_node(conf_formulation, SPC_CF_S_FORMULATION_A_MINERALS, SPC_CF_S_FORMULATION); const YAML::Node& minerals_conf = conf_formulation[SPC_CF_S_FORMULATION_A_MINERALS]; for (auto solid_conf : minerals_conf) { check_mandatory_yaml_node(solid_conf, SPC_CF_S_FORMULATION_A_LABEL, SPC_CF_S_FORMULATION_A_AQUEOUS); check_mandatory_yaml_node(solid_conf, SPC_CF_S_FORMULATION_A_AMOUNT, SPC_CF_S_FORMULATION_A_AQUEOUS); check_mandatory_yaml_node(solid_conf, SPC_CF_S_FORMULATION_A_UNIT, SPC_CF_S_FORMULATION_A_AQUEOUS); const std::string label = solid_conf[SPC_CF_S_FORMULATION_A_LABEL].as(); scalar_t factor = 1.0; units::AmountUnit unit; units::parse_amount_unit(solid_conf[SPC_CF_S_FORMULATION_A_UNIT].as(), unit); // convert to mole concentration is SI if (unit.type == units::AmountUnitType::MoleConcentration or unit.type == units::AmountUnitType::NumberOfMoles) { factor = unit.factor_si; } else if (unit.type == units::AmountUnitType::MassConcentration or unit.type == units::AmountUnitType::Mass ) { const scalar_t molar_mass = molar_mass_solid_unknown_class(label, raw_db); factor = unit.factor_si / molar_mass; } else if (unit.type == units::AmountUnitType::Volume) { const scalar_t molar_volume = molar_volume_solid_unknown_class(label, raw_db); // may raise an error factor = unit.factor_si / molar_volume; } else { throw std::invalid_argument("'"+solid_conf[SPC_CF_S_FORMULATION_A_UNIT].as() + "' is not a valid unit for a solid phase"); } // set to the correct units switch (unit_set.length) { case units::LengthUnit::decimeter: factor *= 1e-3; break; case units::LengthUnit::centimeter: factor *= 1e-6; break; case units::LengthUnit::meter: break; // just to remove warning } const scalar_t mole_conc = factor*solid_conf[SPC_CF_S_FORMULATION_A_AMOUNT].as(); formulation.add_mineral(label, mole_conc); } } // Constraints // =========== // new interface // ------------- void configure_specmicp_constraints( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints ) { // equation for water // ================== // enable conservation of water if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER)) { auto tmp_water = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER); if (tmp_water) { // only one constraints if (conf_constraints.get_optional_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM, false)) { conf_constraints.report_error(YAMLConfigError::InvalidArgument, "Attributes " SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "and " SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM " cannot be set to true at the same time"); } constraints.enable_conservation_water(); } else constraints.disable_conservation_water(); } // or enable saturated system if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) { if (conf_constraints.get_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) { constraints.set_saturated_system(); } } // Aqueous components // ================== // charge keeper // ------------- if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER)) { const auto label = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER); const index_t id = raw_db->get_id_component(label); if (id == no_species) { conf_constraints.report_error( YAMLConfigError::InvalidArgument, "Species " + label + " does not exist in the database." " It can't be the charge keeper." ); } constraints.set_charge_keeper(id); } // Fixed activity // --------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY)) { configure_specmicp_constraints_fixed_activity( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY) ); } // Fixed fugacity // -------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY)) { configure_specmicp_constraints_fixed_fugacity( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY) ); } // Fixed activity // --------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY)) { configure_specmicp_constraints_fixed_molality( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY) ); } // Surface model // ============= if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL)) { const bool is_enabled = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL ); if (is_enabled) { const auto value = conf_constraints.get_required_attribute( SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION); constraints.enable_surface_model(value); } else { constraints.disable_surface_model(); } } } // helpers functions index_t get_id_component( const database::DataContainer* const raw_db, YAMLConfigHandle& conf, const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL ) { const auto label = conf.get_required_attribute(attr_key); const index_t id = raw_db->get_id_component(label); if (id == no_species) { conf.report_error( YAMLConfigError::InvalidArgument, "'" + label + "' is not a valid component." ); } return id; } index_t get_id_gas( const database::DataContainer* const raw_db, YAMLConfigHandle& conf, const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL_GAS ) { const auto label = conf.get_required_attribute(attr_key); const index_t id = raw_db->get_id_gas(label); if (id == no_species) { conf.report_error( YAMLConfigError::InvalidArgument, "'" + label + "' is not a valid gas." ); } return id; } scalar_t get_value(YAMLConfigHandle& conf) { scalar_t value; if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) { value = std::log10(conf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)); } else if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) { value = conf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG); } else { conf.report_error( YAMLConfigError::InvalidArgument, "Expected one of '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '" SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "'." ); } return value; } void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_activity ) { for (uindex_t ind=0; ind(); const index_t id = raw_db->get_id_component(label); if (id == no_species) { throw std::invalid_argument("'"+label+"' is not a valid component (charge_keeper)."); } constraints.set_charge_keeper(id); } // Fixed activity // -------------- if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY]) { configure_specmicp_constraints_fixed_activity( constraints, conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY], raw_db); } // Fixed fugacity // -------------- if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY]) { configure_specmicp_constraints_fixed_fugacity( constraints, conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY], raw_db); } // Fixed molality // -------------- if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY]) { configure_specmicp_constraints_fixed_molality( constraints, conf_constraints[SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY], raw_db); } // Water // ===== if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER]) { if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM] and conf_constraints[SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM].as()) { throw std::invalid_argument("Options '" SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "' and '" SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM "' cannot be used at the same time"); } const bool value = conf_constraints[SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER].as(); if (value) constraints.enable_conservation_water(); else constraints.disable_conservation_water(); } if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM]) { const bool value = conf_constraints[SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM].as(); if (value) constraints.set_saturated_system(); } // Surface model // ============= if (conf_constraints[SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL]) { const bool is_enabled = conf_constraints[SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL].as(); if (is_enabled) { check_mandatory_yaml_node( conf_constraints, SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION, SPC_CF_S_CONSTRAINTS); const scalar_t value = conf_constraints[SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION].as(); constraints.enable_surface_model(value); } else { constraints.disable_surface_model(); } } } void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_activity, const RawDatabasePtr& raw_db ) { for (auto fixact_conf: conf_constraints_fixed_activity) { check_mandatory_yaml_node(fixact_conf, SPC_CF_S_CONSTRAINTS_A_LABEL, SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY); const std::string label = fixact_conf[SPC_CF_S_CONSTRAINTS_A_LABEL].as(); const index_t id = raw_db->get_id_component(label); if (id == no_species) { throw std::invalid_argument("'"+label+"' is not a valid component (fixed_activity)."); } scalar_t value; if (fixact_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT]) { value = std::log10(fixact_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT].as()); } else if (fixact_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG]) { value = fixact_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG].as(); } else { throw std::invalid_argument("Either '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '" SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "' are required in section '" SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY "'."); } constraints.add_fixed_activity_component(id, value); } } void configure_specmicp_constraints_fixed_fugacity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_fugacity, const RawDatabasePtr& raw_db ) { for(auto fixfug_conf: conf_constraints_fixed_fugacity) { check_mandatory_yaml_node(fixfug_conf, SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT, SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY); check_mandatory_yaml_node(fixfug_conf, SPC_CF_S_CONSTRAINTS_A_LABEL_GAS, SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY); const std::string& label_c = fixfug_conf[SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT].as(); const index_t id_c = raw_db->get_id_component(label_c); if (id_c == no_species) { throw std::invalid_argument("'" + label_c + "' is not a valid component (fixed_fugacity)"); } const std::string& label_g = fixfug_conf[SPC_CF_S_CONSTRAINTS_A_LABEL_GAS].as(); const index_t id_g = raw_db->get_id_gas(label_g); if (id_g == no_species) { throw std::invalid_argument("'" + label_g + "' is not a valid gas (fixed_fugacity)"); } scalar_t value; if (fixfug_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT]) { value = std::log10(fixfug_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT].as()); } else if (fixfug_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG]) { value = fixfug_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG].as(); } else { throw std::invalid_argument("Either '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '" SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "' are required in section '" SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY "'."); } constraints.add_fixed_fugacity_gas(id_g, id_c, value); } } void configure_specmicp_constraints_fixed_molality( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_molality, const RawDatabasePtr& raw_db ) { for (auto fixmol_conf: conf_constraints_fixed_molality) { check_mandatory_yaml_node(fixmol_conf, SPC_CF_S_CONSTRAINTS_A_LABEL, SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY); const std::string& label = fixmol_conf[SPC_CF_S_CONSTRAINTS_A_LABEL].as(); const index_t id = raw_db->get_id_component(label); if (id == no_species) { throw std::invalid_argument("'" + label + "' is not a valid component (fixed_molality)"); } scalar_t value; if (fixmol_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT]) { value = std::log10(fixmol_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNT].as()); } else if (fixmol_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG]) { value = fixmol_conf[SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG].as(); } else { throw std::invalid_argument("Either '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '" SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "' are required in section '" SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY "'."); } constraints.add_fixed_molality_component(id, value); } } scalar_t molar_mass_aqueous_unknown_class(const std::string& label, const RawDatabasePtr& raw_db) { scalar_t molar_mass; if (raw_db->get_id_component(label) != no_species) { molar_mass = raw_db->molar_mass_basis(raw_db->get_id_component(label), units::MassUnit::kilogram); } else if (raw_db->get_id_aqueous(label) != no_species) { molar_mass = raw_db->molar_mass_aqueous(raw_db->get_id_aqueous(label), units::MassUnit::kilogram); } else if (raw_db->get_id_compound(label) != no_species) { molar_mass = raw_db->molar_mass_compound(raw_db->get_id_compound(label), units::MassUnit::kilogram); } else { throw std::invalid_argument("'"+label+"is not a valid species"); } return molar_mass; } scalar_t molar_mass_solid_unknown_class(const std::string& label, const RawDatabasePtr& raw_db ) { scalar_t molar_mass; if (raw_db->get_id_mineral(label) != no_species) { molar_mass = raw_db->molar_mass_mineral(raw_db->get_id_mineral(label), units::MassUnit::kilogram); } else if (raw_db->get_id_mineral_kinetic(label) != no_species) { molar_mass = raw_db->molar_mass_mineral_kinetic(raw_db->get_id_mineral_kinetic(label), units::MassUnit::kilogram); } else { throw std::invalid_argument("'"+label+"' is not a valid solid phase."); } return molar_mass; } scalar_t molar_volume_solid_unknown_class(const std::string& label, const RawDatabasePtr& raw_db ) { scalar_t molar_volume = -1; if (raw_db->get_id_mineral(label) != no_species) { molar_volume = raw_db->molar_volume_mineral(raw_db->get_id_mineral(label), units::LengthUnit::meter); } else if (raw_db->get_id_mineral_kinetic(label) != no_species) { molar_volume = raw_db->molar_volume_mineral_kinetic(raw_db->get_id_mineral_kinetic(label), units::LengthUnit::meter); } else { throw std::invalid_argument("'"+label+"' is not a valid solid phase."); } if (molar_volume < 0) { throw std::runtime_error("The molar volume for solid phase '"+label+"' is not defined"); } return molar_volume; } } //end namespace io } //end namespace specmicp diff --git a/src/specmicp_common/io/config_yaml_sections.h b/src/specmicp_common/io/config_yaml_sections.h index 5c751a7..ce0911b 100644 --- a/src/specmicp_common/io/config_yaml_sections.h +++ b/src/specmicp_common/io/config_yaml_sections.h @@ -1,318 +1,319 @@ /* ============================================================================= 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_IO_CONFIGYAMLSECTIONS_H #define SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H //! \brief Define sections in the YAML config files #define SPC_CF_NAME_PATH_ROOT "root" // common options // -------------- #define SPC_CF_A_MAX_ITER "maximum_iterations" #define SPC_CF_A_RES_TOL "residual_tolerance" #define SPC_CF_A_ABS_TOL "absolute_tolerance" #define SPC_CF_A_STEP_TOL "step_tolerance" #define SPC_CF_A_MAX_STEP_LENGTH "maximum_step_length" #define SPC_CF_A_MAX_STEP_MAX_ITER "maximum_step_maximum_iterations" #define SPC_CF_A_FILEPATH "file" #define SPC_CF_A_COMPONENT "component" #define SPC_CF_A_NODES "nodes" #define SPC_CF_A_TYPE "type" #define SPC_CF_V_DEFAULT "default" #define SPC_CF_V_PLUGIN "plugin" #define SPC_CF_A_PLUGIN_FILE "plugin_file" // filepath to a plugin #define SPC_CF_A_PLUGIN_NAME "plugin_name" // name of the plugin to use // user variables // -------------- #define SPC_CF_S_USRVARS "vars" #define SPC_CF_S_USRVARS_SS_INT "int" #define SPC_CF_S_USRVARS_SS_FLOAT "float" // plugin manager // -------------- #define SPC_CF_S_PLUGINS "plugins" #define SPC_CF_S_PLUGINS_A_SEARCHDIRS "dirs" #define SPC_CF_S_PLUGINS_A_TOLOAD "to_load" // Logs // ---- // type of output #define SPC_CF_S_LOGS_V_COUT "cout" #define SPC_CF_S_LOGS_V_CERR "cerr" #define SPC_CF_S_LOGS_V_FILE "file" // the different levels #define SPC_CF_S_LOGS_V_CRITICAL "critical" #define SPC_CF_S_LOGS_V_ERROR "error" #define SPC_CF_S_LOGS_V_WARNING "warning" #define SPC_CF_S_LOGS_V_INFO "info" #define SPC_CF_S_LOGS_V_DEBUG "debug" // runtime logs #define SPC_CF_S_LOGS "logs" #define SPC_CF_S_LOGS_A_LEVEL "level" #define SPC_CF_S_LOGS_A_OUTPUT "output" #define SPC_CF_S_LOGS_A_FILEPATH SPC_CF_A_FILEPATH // configuration logs #define SPC_CF_S_CONF_LOGS "configuration_logs" #define SPC_CF_S_CONF_LOGS_A_OUTPUT "output" #define SPC_CF_S_CONF_LOGS_A_FILE SPC_CF_A_FILEPATH // units // ----- #define SPC_CF_S_UNITS "units" #define SPC_CF_S_UNITS_A_LENGTH "length" #define SPC_CF_S_UNITS_A_MASS "mass" // mesh // ---- #define SPC_CF_S_MESH "mesh" #define SPC_CF_S_MESH_A_TYPE "type" #define SPC_CF_S_MESH_SS_UNIFMESH "uniform_mesh" #define SPC_CF_S_MESH_SS_UNIFMESH_A_DX "dx" #define SPC_CF_S_MESH_SS_UNIFMESH_A_NBNODES "nb_nodes" #define SPC_CF_S_MESH_SS_UNIFMESH_A_SECTION "section" #define SPC_CF_S_MESH_SS_RAMPMESH "ramp_mesh" #define SPC_CF_S_MESH_SS_RAMPMESH_A_MIN_DX "min_dx" #define SPC_CF_S_MESH_SS_RAMPMESH_A_MAX_DX "max_dx" #define SPC_CF_S_MESH_SS_RAMPMESH_A_RAMP_LENGTH "length_ramp" #define SPC_CF_S_MESH_SS_RAMPMESH_A_PLATEAU_LENGTH "length_plateau" #define SPC_CF_S_MESH_SS_RAMPMESH_A_SECTION "section" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH "uniform_axisymmetric" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_RADIUS "radius" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_NBNODES "nb_nodes" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_HEIGHT "height" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_DX "dx" // database // -------- #define SPC_CF_S_DATABASE "database" #define SPC_CF_S_DATABASE_A_PATH "path" #define SPC_CF_S_DATABASE_A_SWAP "swap_components" #define SPC_CF_S_DATABASE_A_SWAP_K_IN "in" #define SPC_CF_S_DATABASE_A_SWAP_K_OUT "out" #define SPC_CF_S_DATABASE_A_REMOVE_GAS "remove_gas" #define SPC_CF_S_DATABASE_A_REMOVE_SOLIDS "remove_solids" #define SPC_CF_S_DATABASE_A_REMOVE_SORBED "remove_sorbeds" #define SPC_CF_S_DATABASE_A_EXTRA_GAS "extra_gas" #define SPC_CF_S_DATABASE_A_EXTRA_SOLIDS "extra_solids" #define SPC_CF_S_DATABASE_A_EXTRA_SORBED "extra_sorbeds" #define SPC_CF_S_DATABASE_A_REMOVE_HALF_CELLS "remove_half_cells" #define SPC_CF_S_DATABASE_A_LIST_SOLIDS_TOKEEP "list_solids_to_keep" // SpecMiCP options // ---------------- #define SPC_CF_S_SPECMICP "specmicp_options" #define SPC_CF_S_SPECMICP_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_SPECMICP_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_SPECMICP_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH #define SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER #define SPC_CF_S_SPECMICP_A_SCALING "enable_scaling" #define SPC_CF_S_SPECMICP_A_NONMONOTONE "enable_nonmonotone_linesearch" #define SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION "factor_descent_direction" #define SPC_CF_S_SPECMICP_A_COND_CHECK "threshold_condition_check" #define SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH "threshold_cycling_linesearch" #define SPC_CF_S_SPECMICP_A_NONIDEAL_TOL "non_ideality_tolerance" #define SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER "non_ideality_"SPC_CF_A_MAX_ITER #define SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC "cutoff_total_concentration" #define SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION "restart_concentration" +#define SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC "restart_water_volume_fraction" #define SPC_CF_S_SPECMICP_A_UNDER_RELAXATION "under_relaxation" // SpecMiCP formulation // -------------------- #define SPC_CF_S_FORMULATION "formulation" #define SPC_CF_S_FORMULATION_A_SOLUTION "solution" #define SPC_CF_S_FORMULATION_A_AQUEOUS "aqueous" #define SPC_CF_S_FORMULATION_A_MINERALS "solid_phases" #define SPC_CF_S_FORMULATION_A_AMOUNT "amount" #define SPC_CF_S_FORMULATION_A_UNIT "unit" #define SPC_CF_S_FORMULATION_A_LABEL "label" // SpecMiCP constraints // -------------------- #define SPC_CF_S_CONSTRAINTS "constraints" #define SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER "charge_keeper" #define SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY "fixed_activity" #define SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY "fixed_molality" #define SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY "fixed_fugacity" #define SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM "saturated_system" #define SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "conservation_water" #define SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL "surface_model" #define SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION "surface_site_concentration" #define SPC_CF_S_CONSTRAINTS_A_LABEL "label" #define SPC_CF_S_CONSTRAINTS_A_AMOUNT "amount" #define SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "amount_log" #define SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT "label_component" #define SPC_CF_S_CONSTRAINTS_A_LABEL_GAS "label_gas" // Variables // --------- #define SPC_CF_S_REACTMICPVARIABLES "variables" #define SPC_CF_S_REACTMICPVARIABLES_A_TYPE SPC_CF_A_TYPE // Boundary conditions // ------------------- #define SPC_CF_S_BOUNDARYCONDITIONS "boundary_conditions" // transport options // ----------------- #define SPC_CF_S_DFPM "transport_options" #define SPC_CF_S_DFPM_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_DFPM_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_DFPM_A_ABS_TOL SPC_CF_A_ABS_TOL #define SPC_CF_S_DFPM_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_DFPM_A_TRSHOLD_STATIONARY "threshold_stationary" #define SPC_CF_S_DFPM_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH #define SPC_CF_S_DFPM_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER #define SPC_CF_S_DFPM_A_SPARSE_SOLVER "sparse_solver" #define SPC_CF_S_DFPM_A_LINESEARCH "linesearch" #define SPC_CF_S_DFPM_A_QUASI_NEWTON "quasi_newton" // ReactMiCP coupling solver // ------------------------- #define SPC_CF_S_REACTMICP "reactmicp_options" #define SPC_CF_S_REACTMICP_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_REACTMICP_A_ABS_TOL SPC_CF_A_ABS_TOL #define SPC_CF_S_REACTMICP_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_REACTMICP_A_GOOD_ENOUGH_TOL "good_enough_tolerance" #define SPC_CF_S_REACTMICP_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_REACTMICP_A_IMPL_UPSCALING "implicit_upscaling" // ReactMiCP output // ---------------- #define SPC_CF_S_REACTOUTPUT "reactmicp_output" #define SPC_CF_S_REACTOUTPUT_A_FILEPATH SPC_CF_A_FILEPATH #define SPC_CF_S_REACTOUTPUT_A_TYPE "type" #define SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5 "hdf5" #define SPC_CF_S_REACTOUTPUT_A_TYPE_V_CSV "csv" #define SPC_CF_S_REACTOUTPUT_A_POROSITY "porosity" #define SPC_CF_S_REACTOUTPUT_A_PH "ph" #define SPC_CF_S_REACTOUTPUT_A_DIFFUSIVITY "diffusivity" #define SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL "volume_fraction_solid" #define SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC "total_aqueous_concentration" #define SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC "total_solid_concentration" #define SPC_CF_S_REACTOUTPUT_A_TOT_CONC "total_concentration" #define SPC_CF_S_REACTOUTPUT_A_SATINDEX "saturation_index" // ReactMiCP simulation stuff // -------------------------- #define SPC_CF_S_SIMULINFO "simulation" #define SPC_CF_S_SIMULINFO_A_NAME "name" #define SPC_CF_S_SIMULINFO_A_OUTPUTPREFIX "output_prefix" #define SPC_CF_S_SIMULINFO_A_PRINTITER "print_iter_info" #define SPC_CF_S_SIMULINFO_A_OUTPUTSTEP "output_step" // ReactMiCP adaptive timestep // --------------------------- #define SPC_CF_S_TIMESTEPPER "timestepper" #define SPC_CF_S_TIMESTEP_A_LOWER_BOUND "minimum_dt" #define SPC_CF_S_TIMESTEP_A_UPPER_BOUND "maximum_dt" #define SPC_CF_S_TIMESTEP_A_RESTART_DT "restart_dt" #define SPC_CF_S_TIMESTEP_A_LOWER_TARGET "lower_iterations_target" #define SPC_CF_S_TIMESTEP_A_UPPER_TARGET "upper_iterations_target" #define SPC_CF_S_TIMESTEP_A_AVERAGE_PARAMETER "average_parameter" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_FAILURE "decrease_factor_if_failure" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_MINIMUM "increase_factor_if_minumum" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_DECREASE "decrease_factor" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_INCREASE "increase_factor" // ReactMiCP Staggers // ------------------ #define SPC_CF_S_STAGGERS "staggers" #define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER "upscaling" #define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE SPC_CF_A_TYPE #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER "chemistry" #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE SPC_CF_A_TYPE #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM "equilibrium" #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS "equilibrium_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER "transport" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION "merge_saturation" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS "merge_aqueous" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS "default_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS "saturation_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS "aqueous_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS "gas_options" // Run ! // ----- #define SPC_CF_S_RUN "run" #define SPC_CF_S_RUN_A_RUNUTNIL "run_until" #endif // SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H