diff --git a/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp b/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp index 5a11aa1..9218910 100644 --- a/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp +++ b/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp @@ -1,354 +1,368 @@ /* ============================================================================= 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_solution_extractor.hpp" #include "specmicp_common/physics/laws.hpp" #include "specmicp_database/database.hpp" namespace specmicp { // a simple cache for the conversion factors struct AdimensionalSystemSolutionExtractor::ExtractorCache { ExtractorCache( const units::UnitsSet& units_set, database::DataContainer* raw_data ): m_scal_factor_mol(units::scaling_factor(units_set.quantity)), m_rho_l(laws::density_water(units_set)), m_scal_molar_vol(raw_data->scaling_molar_volume(units_set)) { const scalar_t scaling_factor_l = units::scaling_factor(units_set.length); m_scal_molar_vol = m_scal_factor_mol / std::pow(scaling_factor_l, 3); m_rt = constants::gas_constant*units::celsius(25.0) * (m_scal_molar_vol); } scalar_t m_scal_factor_mol; scalar_t m_rho_l; scalar_t m_scal_molar_vol; scalar_t m_rt; }; AdimensionalSystemSolutionExtractor::AdimensionalSystemSolutionExtractor( const AdimensionalSystemSolution& solution, RawDatabasePtr raw_data, units::UnitsSet units_set): AdimemsionalSystemNumbering(raw_data), units::UnitBaseClass(units_set), m_solution(solution), m_cache(utils::make_pimpl(units_set, raw_data.get())) {} AdimensionalSystemSolutionExtractor::AdimensionalSystemSolutionExtractor( AdimensionalSystemSolutionExtractor&& other ): AdimemsionalSystemNumbering(other.get_database()), units::UnitBaseClass(std::move(other.get_units())), m_solution(std::move(other.get_solution())), m_cache(std::move(other.m_cache)) { } AdimensionalSystemSolutionExtractor::~AdimensionalSystemSolutionExtractor() = default; scalar_t AdimensionalSystemSolutionExtractor::density_water() const { return m_cache->m_rho_l; } scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_water() const { return density_water()*volume_fraction_water(); } scalar_t AdimensionalSystemSolutionExtractor::pH() const { // find species responsible for pH index_t id = m_data->get_id_component("HO[-]"); if (id != no_species) { return 14+log_activity_component(id); } else { id = m_data->get_id_component("H[+]"); if (id != no_species) return -log_activity_component(id); throw std::runtime_error("No component corresponding to the dissociation of water !"); } } scalar_t AdimensionalSystemSolutionExtractor::total_saturation_minerals() const { return m_solution.main_variables.segment(offset_minerals(), m_data->nb_mineral()).sum(); } scalar_t AdimensionalSystemSolutionExtractor::mole_concentration_mineral(index_t mineral) const { return volume_fraction_mineral(mineral)/( m_data->molar_volume_mineral(mineral)*m_cache->m_scal_molar_vol ); } scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_mineral(index_t mineral) const { return mole_concentration_mineral(mineral)*m_data->molar_mass_mineral(mineral, get_units()); } //! \brief Return the total aqueous concentration scalar_t AdimensionalSystemSolutionExtractor::total_aqueous_concentration(index_t component) const { scalar_t conc = molality_component(component); for (index_t aqueous: m_data->range_aqueous()) { if (m_data->nu_aqueous(aqueous, component) != 0.0) { conc += m_data->nu_aqueous(aqueous, component) * molality_aqueous(aqueous); } } return conc; } //! \brief Return the total solid concentration scalar_t AdimensionalSystemSolutionExtractor::total_solid_concentration(index_t component) const { scalar_t conc = 0; for (index_t mineral: m_data->range_mineral()) { if (m_data->nu_mineral(mineral, component) != 0.0) { conc += m_data->nu_mineral(mineral, component) * mole_concentration_mineral(mineral); } } return conc; } //! \brief Return the total immobile concentration scalar_t AdimensionalSystemSolutionExtractor::total_immobile_concentration(index_t component) const { scalar_t conc = total_solid_concentration(component); const scalar_t conc_w = density_water()*volume_fraction_water(); for (index_t sorbed: m_data->range_sorbed()) { if (m_data->nu_sorbed(sorbed, component) != 0.0) { conc += conc_w * m_data->nu_sorbed(sorbed, component) * molality_sorbed_species(sorbed) / m_cache->m_scal_factor_mol; } } return conc; } scalar_t AdimensionalSystemSolutionExtractor::total_gaseous_concentration(index_t component) const { scalar_t conc = 0; for (index_t gas: m_data->range_gas()) { if (m_data->nu_gas(gas, component) != 0.0) { conc += m_data->nu_gas(gas, component) * fugacity_gas(gas) * 1.01325e5 / m_cache->m_rt; } } return conc; } scalar_t AdimensionalSystemSolutionExtractor::total_concentration( index_t component ) const { return density_water()*volume_fraction_water()*total_aqueous_concentration(component) + total_solid_concentration(component) + volume_fraction_gas_phase()*total_gaseous_concentration(component); } Vector AdimensionalSystemSolutionExtractor::total_concentrations() const { Vector tot_conc(Vector::Zero(m_data->nb_component())); tot_conc(0) = total_concentration(0); for (auto aqc: m_data->range_aqueous_component()) { tot_conc(aqc) = total_concentration(aqc); } return tot_conc; } //! Return the saturation index for 'mineral' scalar_t AdimensionalSystemSolutionExtractor::saturation_index(index_t mineral) const { scalar_t saturation_index = - m_data->logk_mineral(mineral); for (index_t component: m_data->range_aqueous_component()) { if (m_data->nu_mineral(mineral, component) == 0.0) continue; saturation_index += m_data->nu_mineral(mineral, component) * log_activity_component(component); } return saturation_index; } //! Return the saturation index for 'mineral_kinetic' scalar_t AdimensionalSystemSolutionExtractor::saturation_index_kinetic(index_t mineral_kinetic) const { scalar_t saturation_index = - m_data->logk_mineral_kinetic(mineral_kinetic); for (index_t component: m_data->range_aqueous_component()) { if (m_data->nu_mineral_kinetic(mineral_kinetic, component) == 0.0) continue; saturation_index += m_data->nu_mineral_kinetic(mineral_kinetic, component) * log_activity_component(component); } return saturation_index; } scalar_t AdimensionalSystemSolutionExtractor::log_molality_component(index_t component) const { if (component == m_data->water_index()) return std::log10(1.0/m_data->molar_mass_basis(component, get_units())); else if (component == m_data->electron_index()) return -HUGE_VAL; return m_solution.main_variables(dof_component(component)); } scalar_t AdimensionalSystemSolutionExtractor::molality_component(index_t component) const { if (component == m_data->water_index()) return 1.0/m_data->molar_mass_basis(component, get_units()); else if (component == m_data->electron_index()) return 0.0; return pow10(log_molality_component(component))/m_cache->m_scal_factor_mol; } scalar_t AdimensionalSystemSolutionExtractor::molality_aqueous(index_t aqueous) const { return m_solution.secondary_molalities(dof_aqueous(aqueous))/m_cache->m_scal_factor_mol; } scalar_t AdimensionalSystemSolutionExtractor::total_mass_concentration() const { return (mass_concentration_water() + total_solid_mass_concentration()); } scalar_t AdimensionalSystemSolutionExtractor::total_solid_mass_concentration() const { scalar_t mass = 0.0; for (index_t mineral: m_data->range_mineral()) { mass += mass_concentration_mineral(mineral); } return mass; } +scalar_t AdimensionalSystemSolutionExtractor::charge_balance() const +{ + scalar_t balance = 0; + + for (index_t i: m_data->range_aqueous_component()) { + balance += m_data->charge_component(i)*molality_component(i); + } + for (index_t j: m_data->range_aqueous()) { + balance += m_data->charge_aqueous(j)*molality_aqueous(j); + } + + return balance; +} + // ########################### // // // // Modificator // // // // ########################### // void AdimensionalSystemSolutionModificator::scale_total_concentration( index_t component, scalar_t new_value) { const scalar_t old_value = total_solid_concentration(component); const scalar_t factor = new_value/old_value; m_nonconst_solution.main_variables.segment(dof_mineral(0), m_data->nb_mineral()) *= factor; } void AdimensionalSystemSolutionModificator::remove_solids() { m_nonconst_solution.main_variables.segment(dof_mineral(0), m_data->nb_mineral()).setZero(); } Vector AdimensionalSystemSolutionModificator::set_minerals_kinetics(std::vector& list_species) { index_t nb_kinetics = list_species.size(); index_t nb_new_mineral = m_data->nb_mineral() - nb_kinetics; std::vector minerals_to_keep; minerals_to_keep.reserve(nb_new_mineral); std::vector new_kinetics_index(nb_kinetics, no_species); Vector saturation_kinetics(nb_kinetics); index_t new_ind_eq = offset_minerals(); index_t new_ind_kin = 0; index_t tot_ind_kin = m_data->nb_mineral_kinetic(); // ###TODO optimize for (index_t mineral: m_data->range_mineral()) { auto is_kin = std::find(list_species.begin(), list_species.end(), mineral); // If mineral is still at equilibrium if (is_kin == list_species.end()) { minerals_to_keep.push_back(mineral); m_nonconst_solution.main_variables(new_ind_eq) = volume_fraction_mineral(mineral); ++new_ind_eq; } // If mineral is governed by kinetics else { saturation_kinetics(new_ind_kin) = volume_fraction_mineral(mineral); ++new_ind_kin; // save the new index (index in the kinetics vector) // The order is conserved ! new_kinetics_index[is_kin - list_species.begin()] = tot_ind_kin; ++tot_ind_kin; } } // change the database database::Database dbhandler(m_data); dbhandler.minerals_keep_only(minerals_to_keep); specmicp_assert(new_ind_eq == total_dofs()); // update the list of species list_species.swap(new_kinetics_index); // resize m_nonconst_solution.main_variables.conservativeResize(total_dofs()); return saturation_kinetics; } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp index 7ea051f..aa97cb1 100644 --- a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp +++ b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp @@ -1,353 +1,356 @@ /* ============================================================================= 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_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP #include "adimensional_system_solution.hpp" #include "adimensional_system_numbering.hpp" #include "specmicp_common/types.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/physics/constants.hpp" #include "specmicp_common/pimpl_ptr.hpp" //! \file adimensional_system_solution_extractor.hpp Obtain data from AdimensionalSystemSolution namespace specmicp { //! \class AdimensionalSystemSolutionExtractor //! \brief Obtain physical variables from AdimensionalSystemSolution //! //! This class can't modify the solution, use specmicp::AdimensionalSystemSolutionModificator for this. //! //! \ingroup specmicp_api //! \sa specmicp::AdimensionalSystemSolutionModificator class SPECMICP_DLL_PUBLIC AdimensionalSystemSolutionExtractor: public AdimemsionalSystemNumbering, public units::UnitBaseClass { public: // public members // ############### //! \brief Constructor //! //! \param solution Reference to a solution from the AdimensionalSystemSolution //! \param thedatabase shared_ptr to the database //! \param units_set Units used by the solver //! //! \warning This is just a reader, it does not store the information AdimensionalSystemSolutionExtractor( const AdimensionalSystemSolution& solution, RawDatabasePtr thedatabase, units::UnitsSet units_set); //! \brief Move constructor AdimensionalSystemSolutionExtractor(AdimensionalSystemSolutionExtractor&& other); //! \brief Destructor ~AdimensionalSystemSolutionExtractor(); // Primary variables // ================= // water // ------ //! \brief Return the volume fraction of water scalar_t volume_fraction_water() const { return m_solution.main_variables(dof_water());} //! \brief Return the density of water scalar_t density_water() const; // electron // --------- //! \brief Return the log of the activity electron scalar_t log_activity_electron() const { return m_solution.main_variables(dof_electron()); } //! \brief Return the activity of the electron scalar_t activity_electron() const { return pow10(log_activity_electron()); } //! \brief Return the pE scalar_t pE() const { return -log_activity_electron(); } //! \brief The electric potential scalar_t Eh() const { return std::log(10.0)*constants::gas_constant*units::celsius(25.0)/constants::faraday_constant*pE(); } // component // --------- //! \brief Return the log_10 of the molality of 'component' scalar_t log_molality_component(index_t component) const; //! \brief Return the molality of 'component' scalar_t molality_component(index_t component) const; // solid phases // ------------ //! \brief Return the volume fraction of 'mineral' scalar_t volume_fraction_mineral(index_t mineral) const { return m_solution.main_variables(dof_mineral(mineral)); } // surface // ------- //! \brief Return the concentration of free surface scalar_t log_free_surface_concentration() const { return m_solution.main_variables(dof_surface()); } //! \brief Return the concentration of free surface scalar_t free_surface_concentration() const { return pow10(m_solution.main_variables(dof_surface())); } // Secondary variables // =================== //! \brief Return the ionic strength of the solution scalar_t ionic_strength() const {return m_solution.ionic_strength;} // component // --------- //! \brief Return the log_10 of the activity coefficient of 'component' scalar_t log_activity_coefficient_component(index_t component) const { return m_solution.log_gamma(dof_component_gamma(component)); } //! \brief Return the activity coefficient of 'component' scalar_t activity_coefficient_component(index_t component) const { return pow10(log_activity_coefficient_component(component)); } //! \brief Return the log_10 of the activity of 'component scalar_t log_activity_component(index_t component) const { return log_molality_component(component) + log_activity_coefficient_component(component); } //! \brief Return the activity of 'component' scalar_t activity_component(index_t component) const { return pow10(log_activity_component(component)); } // aqueous species // --------------- //! \brief Return the molality of secondary specis 'aqueous' scalar_t molality_aqueous(index_t aqueous) const; //! \brief Return the log10 of the activity ocefficient of secondary species 'aqueous' scalar_t log_activity_coefficient_aqueous(index_t aqueous) const { return m_solution.log_gamma(dof_aqueous_gamma(aqueous)); } //! \brief Return the activity coefficient of secondary species 'aqueous' scalar_t activity_coefficient_aqueous(index_t aqueous) const { return pow10(log_activity_coefficient_aqueous(aqueous)); } //! \brief Return the activity of secondary species 'aqueous' scalar_t activity_aqueous(index_t aqueous) const { return activity_coefficient_aqueous(aqueous)*molality_aqueous(aqueous); } // gas fugacity // ------------- //! \brief Return fugacity for 'gas' scalar_t fugacity_gas(index_t gas) const { return m_solution.gas_fugacities(gas); } // Sorbed species // -------------- //! \brief Return sorbed species molalities for 'sorbed' scalar_t molality_sorbed_species(index_t sorbed) const { return m_solution.sorbed_molalities(sorbed); } // Tertiary variables // ================== // Water // ----- //! \brief Return the saturation of water scalar_t saturation_water() const {return volume_fraction_water()/porosity();} //! \brief Return the mass concentration of water scalar_t mass_concentration_water() const; //! \brief Return the pH of the solution scalar_t pH() const; // Component // --------- //! \brief Return the total aqueous concentration w.r.t solution (x*mol/kg) scalar_t total_aqueous_concentration(index_t component) const; //! \brief Return the total solid concentration w.r.t to REV(x*mol/m^3) scalar_t total_solid_concentration(index_t component) const; //! \brief Return the total immobile concentration w.r.t ro REV (x*mol/m^3) scalar_t total_immobile_concentration(index_t component) const; //! \brief Return the total gaseous concentration w.r.t to gas phase (x*mol/m^3) scalar_t total_gaseous_concentration(index_t component) const; //! \brief Return the total concentration vector (x*mol/m^3) scalar_t total_concentration(index_t component) const; //! \brief Return the vector of total concentration (x*mol/m^3) Vector total_concentrations() const; // Gas phase // --------- //! \brief Return the saturation of the gas phase scalar_t saturation_gas_phase() const {return 1-saturation_water();} //! \brief Return the volume fraction occupied by the gas phase scalar_t volume_fraction_gas_phase() const { return porosity() - volume_fraction_water(); } // Component // --------- // Solid phases // ------------ //! \brief Return the concentration of 'mineral' //! //! \param mineral Index of the mineral (in the database) scalar_t mole_concentration_mineral(index_t mineral) const; //! \brief Return the mass concentration of 'mineral' //! //! \param mineral Index of the mineral (in the database) scalar_t mass_concentration_mineral(index_t mineral) const; //! \brief Return the total saturation of all minerals //! Deprecated, use the volume fraction function instead scalar_t total_saturation_minerals() const; //! \brief Return the volume fraction occupied by the minerals scalar_t volume_fraction_minerals() const {return total_saturation_minerals();} //! \brief Return the inert volume fraction scalar_t volume_fraction_inert() const {return m_solution.inert_volume_fraction;} //! \brief Return the porosity scalar_t porosity() const {return 1 - total_saturation_minerals() - volume_fraction_inert();} //! \brief Return the saturation index for 'mineral' //! //! \param mineral Index of the mineral (in the database) scalar_t saturation_index(index_t mineral) const; //! \brief Return the saturation index for 'mineral_kinetic' //! //! \param mineral_kinetic Index of the mineral governed by the kinetic (in the database) scalar_t saturation_index_kinetic(index_t mineral_kinetic) const; //! \brief Return the total mass concentration of the system (liquid + solid at equilibrium) scalar_t total_mass_concentration() const; //! \brief Return the total mass concentration of the solids at equilibrium) scalar_t total_solid_mass_concentration() const; + //! \brief Return the charge molality + scalar_t charge_balance() const; + //! \brief Return the vector of main variables, can be used to initialize the solver Vector get_main_variables() const {return m_solution.main_variables;} //! \brief Return a copy of the solution AdimensionalSystemSolution get_solution() const {return m_solution;} protected: // private attributes // ################## const AdimensionalSystemSolution& m_solution; //!< Reference to the solution //! \brief A cache to store some computations struct SPECMICP_DLL_LOCAL ExtractorCache; utils::pimpl_ptr m_cache; //!< Cache to avoid recomputing stuff }; //! \class AdimensionalSystemSolutionModificator //! \brief Special class to modify (scale) the solution //! //! The set of modification are limited to straightforward modification that does not modify the equilibrium //! //! \sa AdimensionalSystemExtractor //! \ingroup specmicp_api class SPECMICP_DLL_PUBLIC AdimensionalSystemSolutionModificator: public AdimensionalSystemSolutionExtractor { public: //! \brief Constructor AdimensionalSystemSolutionModificator( AdimensionalSystemSolution& solution, RawDatabasePtr thedatabase, units::UnitsSet units_set): AdimensionalSystemSolutionExtractor(solution, thedatabase, units_set), m_nonconst_solution(solution) {} //! \brief Scale the solid phases with respect to the total solid concentration of a component //! //! \param component Index of the component (in the database) //! \param new_value New value of the total solid concentration void scale_total_concentration(index_t component, scalar_t new_value); //! \brief Remove the solid phases void remove_solids(); //! \brief Set some species to be considered as governed by kinetics //! //! \param[in,out] list_species list of minerals to flag as kinetics, //! The list will contain the new indexes of the corresponding solid phases; //! \return the vector of saturation //! //! This function modify the databse Vector set_minerals_kinetics(std::vector& list_species); private: AdimensionalSystemSolution& m_nonconst_solution; //!< Brief Reference to the solution }; } // end namespace specmicp #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP