diff --git a/data/cemdata_specmicp.js b/data/cemdata_specmicp.js index 0dc00d3..7799d57 100644 --- a/data/cemdata_specmicp.js +++ b/data/cemdata_specmicp.js @@ -1,515 +1,515 @@ // This is a manual and partial translation of the cemdataO7 database // The list of references is at the end of the database // The database is in the JSON format, to be read easily by both human and machine.z { "Metadata": { - "version": "0.0.2", - "last-modification": "04/30/2014" + "version": "0.0.3", + "last-modification": "06/05/2014" }, "Basis": [ { "label": "H2O", "molar_mass": 18.015 // g/mol }, { "label": "H[+]", "molar_mass": 1.008, "activity": { "a": 9.00, "b": 0.00 } }, { "label": "Al[3+]", "molar_mass": 26.982, "activity": { "a": 6.65, "b": 0.19 } }, { "label": "HCO3[-]", "molar_mass": 61.016, "activity": { "a": 5.40, "b": 0.00 } }, { "label": "Ca[2+]", "molar_mass": 40.078, "activity": { "a": 4.86, "b": 0.15 } }, { "label": "Cl[-]", "molar_mass": 35.453, "activity": { "a": 3.71, "b": 0.01 } }, { "label": "Na[+]", "molar_mass": 22.99, "activity": { "a": 4.32, "b": 0.06 } }, { "label": "Si(OH)4", "molar_mass": 96.114, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "SO4[2-]", "molar_mass": 96.063, "activity": { "a": 5.31, "b": -0.07 } } ], "Aqueous": [ { "label": "HO[-]", "composition": "H2O, - H[+]", "log_k": 13.9995, "activity": { "a": 10.65, "b": 0.00 } }, { "label": "AlOH[2+]", "composition": "Al[3+], H2O, -H[+]", "log_k": 4.9573, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "Al(OH)2[+]", "composition": "Al[3+], 2 H2O, -2 H[+]", "log_k": 10.5943, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "Al(OH)3", "composition": "Al[3+], 3 H2O, -3 H[+]", "log_k": 16.4328, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "Al(OH)4[-]", "composition": "Al[3+], 4 H2O, -4 H[+]", "log_k": 22.8797, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "AlSO4[+]", "composition": "Al[3+], SO4[2-]", "log_k": -3.900, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "Al(SO4)2[-]", "composition": "Al[3+], 2 SO4[2-]", "log_k": -5.9000, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "HSO4[-]", "composition": "SO4[2-], H[+]", "log_k": -1.9878, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CO2", "composition": "-H2O, H[+], +HCO3[-]", "log_k": -6.3519, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "CO3[2-]", "composition": "HCO3[-],- H[+]", "log_k": 10.3289, "activity": { "a": 5.40, "b": 0.00 } }, { "label": "CaSO4", "composition": "Ca[2+], SO4[2-]", "log_k": -2.300, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "SiO(OH)3[-]", "composition": "Si(OH)4, -H[+]", "log_k": 9.8100, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "SiO2(OH)2[2-]", "composition": "Si(OH)4, -2 H[+]", "log_k": 23.1400, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "AlSiO(OH)3[2+]", "composition": "Al[3+], SiO(OH)3[-]", "log_k": -7.400, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "Al(OH)6SiO[-]", "composition": "Al(OH)4[-], Si(OH)4, -H2O", "log_k": -3.600, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaCO3", "composition": "HCO3[-], Ca[2+], -H[+]", "log_k": 7.1048, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "CaHCO3[+]", "composition": "HCO3[-], Ca[2+]", "log_k": -1.1057, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaOH[+]", "composition": "Ca[2+], H2O, -H[+]", "log_k": 12.7800, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaSiO(OH)3[+]", "composition": "Ca[2+], SiO(OH)3[-]", "log_k": -1.200, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaSiO2(OH)2", "composition": "Ca[2+], SiO2(OH)2[2-]", "log_k": -4.60, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "NaCO3[-]", "composition": "Na[+], -H[+], HCO3[-]", "log_k": 9.0590, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "NaHCO3", "composition": "Na[+], HCO3[-]", "log_k": 0.2500, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "NaOH", "composition": "Na[+], H2O, -H[+]", "log_k": 14.18, "activity": { "a": 0.00, "b": 0.10 } }, { "label": "NaSO4[-]", "composition": "Na[+], SO4[2-]", "log_k": -0.7, "activity:": { "a": 4.00, "b": 0.00 } }, { "label": "HSO4[-]", "composition": "H[+], SO4[2-]", "log_k": -1.9878, "activity": { "a": 4.00, "b": 0.00 } } ], "Minerals": [ // Hydrates { "label": "Portlandite", "composition": "Ca[2+], +2 H2O, -2 H[+]", "log_k": 22.800, "molar_volume": 33, "references": ["Hummel2002"] }, { "label": "CSH,jennite", "composition": "1.6667 Ca[2+], 2.3334 HO[-], SiO(OH)3[-], -0.5567 H2O", "log_k": -13.17, "molar_volume": 78, "references": ["Lothenbach2008"] }, { "label": "CSH,tobermorite", "composition": "0.8333 Ca[2+], 0.6666 HO[-], SiO(OH)3[-], -0.5 H2O", "log_k": -8.00, "molar_volume": 59, "references": ["Lothenbach2008"] }, { "label": "Al(OH)3,am", "composition": "Al(OH)4[-], -HO[-]", "log_k": 0.24, "molar_volume": 32, "references": ["Lothenbach2008"] }, { "label": "Tricarboaluminate", "composition": "6 Ca[2+], 2 Al(OH)4[-], 3 CO3[2-], 4 HO[-], 26 H2O", "log_k": -46.5, "molar_volume": 650, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "C3AH6", "composition": "3 Ca[2+], 2 Al(OH)4[-], 4 HO[-]", "log_k": -20.84, "molar_volume": 150, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "C4AH13", "composition": "4 Ca[2+], 2 Al(OH)4[-], 6 HO[-], 6 H2O", "log_k": -25.40, "molar_volume": 274, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "C2AH8", "composition": "2 Ca[2+], 2 Al(OH)4[-], 2 HO[-], 3 H2O", "log_k": -13.56, "molar_volume": 184, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Monosulfoaluminate", "composition": "4 Ca[2+], 2 Al(OH)4[-], SO4[2-], 4 HO[-], 6 H2O", "log_k": -29.26, "molar_volume": 309, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Monocarboaluminate", "composition": "4 Ca[2+], 2 Al(OH)4[-], CO3[2-], 4 HO[-], 5 H2O", "log_k": -31.47, "molar_volume": 262, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Hemicarboaluminate", "composition": "4 Ca[2+], 2 Al(OH)4[-], 0.5 CO3[2-], 5 HO[-], 5.5 H2O", "log_k": -29.13, "molar_volume": 285, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Straetlingite", "composition": "2 Ca[2+], 2 Al(OH)4[-], SiO(OH)3[-], HO[-], 2 H2O", "log_k": -19.70, "molar_volume": 216, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "CAH10", "composition": "Ca[2+], 2 Al(OH)4[-], 6 H2O", "log_k": -7.50, "molar_volume": 194, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Gypsum", "composition": "Ca[2+], SO4[2-], 2 H2O", "log_k": -4.5809, "molar_volume": 74, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Ettringite", "composition": "6 Ca[2+], 2 Al(OH)4[-], 3 SO4[2-], 4 HO[-], 26 H2O", "log_k": -44.9, "molar_volume": 707, "references": ["Lothenbach2008","Matschei2007"] }, { "label": "Thaumasite", "composition": "6 Ca[2+], 2 SiO(OH)3[-], 2 SO4[2-], 2 CO3[2-], 2 HO[-], 26 H2O", "log_k": -49.4, "molar_volume": 663, "references": ["Schmidt2008"] }, // Anhydrous { "label": "Calcite", "composition": "Ca[2+], HCO3[-], -H[+]", "log_k": 1.8490, "molar_volume": 37, "references": ["Hummel2002"] }, { "label": "SiO2,am", "composition": "SiO(OH)3[-], - H2O, - HO[-]", "log_k": 1.476, "molar_volume": 29, "references": ["Lothenbach2008"] }, // Others { "label": "Kaolinite", "composition": "2 Al[3+], 2 Si(OH)4, H2O, -6 H[+]", "log_k": 7.4350, "flag_kinetic": "true", "references": ["Hummel2002"] }, { "label": "Gibbsite", "composition": "Al[3+], 3 H2O, -3 H[+]", "log_k": 7.7561, "flag_kinetic": "true", "references": ["Hummel2002"] }, { "label": "C3S", "composition": "3 Ca[2+], Si(OH)4, 6 HO[-]", "log_k": -22.0, "flag_kinetic": "true", "molar_volume": 73, "references": ["Nicoleau2013"] }, { "label": "C2S", "composition": "2 Ca[2+], Si(OH)4, 4HO[-]", "log_k": -17.4, "molar_volume": 52, "flag_kinetic": ["Nicoleau2013"] }, { "label": "C3A", "composition": "3 Ca[2+], 8 HO[-], -6 H2O, 2 Al(OH)2[+]", "log_k": 50, "molar_volume": 89, "flag_kinetic": "true" }, { "label": "Anhydrite", "composition": "Ca[2+], SO4[2-]", "log_k": -4.3575, "flag_kinetic": "true", "molar_volume": 46, "references": ["Lothenbach2008","Matschei2007"] } ], "References" : [ { "key": "Hummel2002", "description": "Hummel W., Berner U., Curti E., Pearson F.J. & Thoenen T. (2002): Nagra/PSI Chemical Thermodynamic Data Base 01/01. Universal Publishers/uPUBLISH.com USA, available from: http://www.upublish.com/books/hummel.htm. Also issued as Nagra Technical Report NTB 02-16, Nagra, Wettingen, Switzerland." }, { "key": "Pearson1991", "description": "Pearson F.J. & Berner U. (1991): Nagra Thermochemical Data Base I. Core Data. Nagra Technical Report NTB 91-17, Nagra, Wettingen, Switzerland." }, { "key": "Pearson1992", "description": " Pearson F.J., Berner U. & Hummel W. (1992): Nagra Thermochemical Data Base II. Supplemental Data 05/92. Nagra Technical Report NTB 91-18, Nagra, Wettingen, Switzerland." }, { "key": "Matschei2007", "description": "Matschei, Thomas, Lothenbach, Barbara, Glasser, Fredrik P.: Thermodynamic properties of Portland cement hydrates in the system CaO–Al2O3–SiO2–CaSO4–CaCO3–H2O, Cement and Concrete Research 37(10), 1379–1410, 2007" }, { "key": "Lothenbach2008", "description": "Lothenbach, Barbara, Matschei, Thomas, Möschner, Göril, Glasser, Fred P.: Thermodynamic modelling of the effect of temperature on the hydration and porosity of Portland cement, Cement and Concrete Research 38(1), 1–18, 2008" }, { "key": "Schmidt2008", "description": "Schmidt, Thomas, Lothenbach, Barbara, Romer, Michael, Scrivener, Karen, Rentsch, Daniel, Figi, Renato: A thermodynamic and experimental study of the conditions of thaumasite formation, Cement and Concrete Research 38(3), 337–349, 2008" }, { "key": "Nicoleau2013", "description": "Nicoleau, L., Nonat, A., Perrey, D.: The di- and tricalcium silicate dissolutions, Cement and Concrete Research 47(0), 14–30, 2013" } ] } diff --git a/data/specmicp_database.js b/data/specmicp_database.js index 4b8fb48..0175dd1 100644 --- a/data/specmicp_database.js +++ b/data/specmicp_database.js @@ -1,202 +1,211 @@ { "Metadata": { "version": "0.0.1", - "last-modification": "01/12/2013" + "last-modification": "06/05/2014" }, "Basis": [ { - "label": "H2O" + "label": "H2O", + "molar_mass": 18.015 // g/mol + }, { "label": "HO[-]", + "molar_mass": 1.008, "activity": { "a": 10.65, "b": 0 } }, { "label": "Ca[2+]", + "molar_mass": 40.078, "activity": { "a": 4.86, "b": 0.15 } }, { "label": "SiO(HO)3[-]", + "molar_mass": 95.008, "activity": { "a": 10.65, "b": 0.00 } }, { "label": "SO4[2-]", + "molar_mass": 96.063, "activity": { "a": 5.31, "b": -0.07 } }, { "label": "HCO3[-]", + "molar_mass": 61.016, "activity": { "a": 5.40, "b": 0.00 } }, { "label": "Na[+]", + "molar_mass": 22.99, "activity": { "a": 4.32, "b": 0.06 } }, { "label": "Cl[+]", + "molar_mass": 35.453, "activity": { "a": 3.71, "b": 0.01 } } ], "Aqueous": [ { "label": "H[+]", "composition": "H2O, - HO[-]", "log_k": 13.9995, "activity": { "a": 9.00, "b": 0.00 } }, { "label": "Si(OH)4", "composition": "SiO(HO)3[-], H2O, - HO[-]", "log_k": 4.8595 }, { "label": "SiO2(OH)2[2-]", - "composition": "SiO(HO)3[-], H2O, - HO[-]", + "composition": "SiO(HO)3[-], - H2O, + HO[-]", "log_k": 0.1005, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaSO4(aq)", "composition": "Ca[2+], SO4[2-]", "log_k": -2.3 }, { - "label": "CaOH+", + "label": "CaOH[+]", "composition": "Ca[2+], HO[-]", "log_k": -1.2195, "activity":{ "a": 4.00, "b": 0.00 } }, { "label": "CO3[2-]", "composition": "HCO3[-], - H2O, HO[-]", "log_k": -3.6706, "activity": { "a": 5.40, "b": 0.00 } }, { "label": "CO2(aq)", "composition": "HCO3[-], - HO[-]", "log_k": 7.6476 }, { "label": "CaCO3(aq)", "composition": "Ca[2+], HCO3[-], HO[-], -H2O", "log_k": -6.8947 }, { "label": "CaHCO3[+]", "composition": "Ca[2+], HCO3[-]", "log_k": -1.1057, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaSiO(HO)3[+]", "composition": "Ca[2+], SiO(HO)3[-]", "log_k": -1.2, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "CaSiO2(HO)2(aq)", "composition": "Ca[2+], SiO(HO)3[-], HO[-], - H2O", "log_k": -4.4995 }, { "label": "NaOH(aq)", "composition": "Na[+], HO[-]", "log_k": 0.1805 }, { "label": "NaSO4[-]", "composition": "Na[+], SO4[2-]", "log_k": -0.7, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "NaCO3[-]", "composition": "Na[+], HO[-], HCO3[-], - H2O", "log_k": -4.9405, "activity": { "a": 4.00, "b": 0.00 } }, { "label": "NaHCO3(aq)", "composition": "Na[+], HCO3[-]", "log_k": 0.25 } ], "Minerals": [ { "label": "Portlandite", "composition": "Ca[2+], 2 HO[-]", "log_k": -5.1995 }, { "label": "SiO2,am", "composition": "SiO(HO)3[-], - H2O, - HO[-]", "log_k": 1.476 }, { "label": "CSH,jennite", - "composition": "1.6667 Ca[2+], 2.3333 HO[-], SiO(HO)3[-], -0.5567 H2O", + "composition": "1.6667 Ca[2+], 2.3334 HO[-], SiO(HO)3[-], -0.5567 H2O", "log_k": -13.17 }, { "label": "CSH,tobermorite", - "composition": "0.8333 Ca[2+], 0.6667 HO[-], SiO(HO)3[-], -0.5 H2O", + "composition": "0.8333 Ca[2+], 0.6666 HO[-], SiO(HO)3[-], -0.5 H2O", "log_k": -8.00 }, { "label": "Calcite", "composition": "Ca[2+], HCO3[-], HO[-], -H2O", "log_k": -12.1505 }, { "label": "Gypsum", "composition": "Ca[2+], SO4[2-], 2 H2O", "log_k": -4.5809 } ] } diff --git a/src/odeint/odeint_types.hpp b/src/odeint/odeint_types.hpp index 238c349..7cc0bd1 100644 --- a/src/odeint/odeint_types.hpp +++ b/src/odeint/odeint_types.hpp @@ -1,67 +1,77 @@ /*------------------------------------------------------- - Module : odeint - File : odeint_types.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_ODEINT_ODEINTTYPES_HPP #define SPECMICP_ODEINT_ODEINTTYPES_HPP //! \file odeint_types.hpp common types for odeint #include #include namespace specmicp { namespace odeint { //! A scalar using scalar_t = double; namespace internal { -//! implementation for the vector type (which may not be a vector for dim=1) +//! \brief implementation for the vector type (which may not be a vector for dim=1) template struct vector_trait { - using type = Eigen::Matrix; + using type = Eigen::Matrix; }; +//! Specialisation for dimension 1 (scalar) template <> struct vector_trait<1> { using type = scalar_t; }; +//! Specialization for dimension 2 (fixed-size vector) +template <> +struct vector_trait<2> +{ + using type = Eigen::Matrix; +}; + +//! Specialization for dimension 3 (fixed-size vector) +template <> +struct vector_trait<3> +{ + using type = Eigen::Matrix; +}; + + } // end namespace internal //! The vector type template using vector_t = typename internal::vector_trait::type; -//! +//! Return the maximum coefficient of the vector template scalar_t max_coeff(const vector_t& vec) {return vec.maxCoeff();} template <> scalar_t max_coeff<1>(const vector_t<1>& vec) {return vec;} -template -void copy(vector_t& receiver, const vector_t& donor) -{ - receiver = donor; -} - //! Type of the function f in dy/dx=f(x,y) template using rhs_f = std::function&, vector_t&)>; } // end namespace odeint } // end namespace specmicp #endif // SPECMICP_ODEINT_ODEINTTYPES_HPP diff --git a/src/specmicp/equilibrium_data.cpp b/src/specmicp/equilibrium_data.cpp index a7cc123..bb293c6 100644 --- a/src/specmicp/equilibrium_data.cpp +++ b/src/specmicp/equilibrium_data.cpp @@ -1,70 +1,101 @@ /*------------------------------------------------------- - Module : specmicp - File : equilibrium_data.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "equilibrium_data.hpp" #include "database/module.hpp" namespace specmicp { //! return log(IAP) for mineral m double EquilibriumState::logIAP(int m) { double logiap = 0.0; for (int i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral(m, i); return logiap; } +double EquilibriumState::logIAP(const std::string& label_mineral) +{ + database::DatabaseModule namer(m_data); + int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); + return logIAP(idm); +} + + +double EquilibriumState::saturation_index(const std::string& label_mineral) +{ + database::DatabaseModule namer(m_data); + int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); + return logIAP(idm) - m_data->logk_mineral(idm); +} + //! return log(IAP) for mineral m double EquilibriumState::logIAP_kinetic(int m) { double logiap = 0.0; for (int i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral_kinetic(m, i); return logiap; } +double EquilibriumState::logIAP_kinetic(const std::string& label_mineral) +{ + database::DatabaseModule namer(m_data); + int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals_kinetic); + return logIAP_kinetic(idm); +} + + +double EquilibriumState::saturation_index_kinetic(const std::string& label_mineral) +{ + database::DatabaseModule namer(m_data); + int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); + return logIAP_kinetic(idm) - m_data->logk_mineral_kinetic(idm); +} + //! Return total concentration of aqueous species for component i double EquilibriumState::total_aqueous_concentration_component(int i) { double totconc = POW10(m_main_variables(i)); for (int j=0; jnb_aqueous; ++j) { if (m_data->nu_aqueous(j,i) != 0) totconc += m_data->nu_aqueous(j,i)*m_concentration_secondary(j); } return totconc; } void EquilibriumState::total_aqueous_concentrations(Eigen::VectorXd& totconc) { assert(totconc.rows() == m_data->nb_component); totconc(0) = m_main_variables(0); for (int i=1; inb_component; ++i) totconc(i) = total_aqueous_concentration_component(i); } double EquilibriumState::pH() { + // search the information about "H[+]" database::DatabaseModule namer(m_data); double id = namer.component_label_to_id("H[+]"); - if (id != database::DatabaseModule::no_species) return m_main_variables(id); + if (id != database::DatabaseModule::no_species) return -m_main_variables(id); else { id = namer.aqueous_label_to_id("H[+]"); assert(id != database::DatabaseModule::no_species); return -std::log10(m_concentration_secondary(id)); } } double EquilibriumState::mass_mineral(int m) { return m_data->molar_mass_mineral(m)*m_main_variables(m_data->nb_component+m); } } // end namespace specmicp diff --git a/src/specmicp/equilibrium_data.hpp b/src/specmicp/equilibrium_data.hpp index f8c6687..0297a3d 100644 --- a/src/specmicp/equilibrium_data.hpp +++ b/src/specmicp/equilibrium_data.hpp @@ -1,107 +1,115 @@ /*------------------------------------------------------- - Module : specmicp - File : equilibrium_data.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP #define SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP //! \file equilibrium_data.hpp store results from equilibrium computation #include "database/database.hpp" #define POW10(m) std::pow(10.0, m) namespace specmicp { class EquilibriumState { EIGEN_MAKE_ALIGNED_OPERATOR_NEW public: EquilibriumState() {} EquilibriumState(const Eigen::VectorXd& main_vars, const Eigen::VectorXd& secondary, const Eigen::VectorXd& loggamma, double ionic_strength, std::shared_ptr data): m_main_variables(main_vars), m_concentration_secondary(secondary), m_ionic_strength(ionic_strength), m_loggamma(loggamma), m_data(data) {} // Main variables // ============== //! Return the mass of free water in the system double mass_water() {return m_main_variables(0);} //! Return the molality of component i double molality_component(int i) {return POW10(m_main_variables(i));} //! Return the number of moles of the mineral at equilibrium 'm' double moles_mineral(int m) {return m_main_variables(m_data->nb_component+m);} //! Return the mass of a mineral (kg) double mass_mineral(int m); // molality // ======== // return the pH double pH(); //! Return the molality of secondary species j double molality_secondary(int j) {return m_concentration_secondary(j);} //! Return total concentration of aqueous species for component i double total_aqueous_concentration_component(int i); //! Compute total aqeuous concentration for each component void total_aqueous_concentrations(Eigen::VectorXd& totconc); // Activity // ======== //! Return activity coefficient of component i double activity_coefficient_component(int i) { assert(i > 0); return POW10(m_loggamma(i)); } //! Return activity coefficient of secondary species j double activity_coefficient_secondary(int j) {return POW10(m_loggamma(j+m_data->nb_component));} //! Return activity of component i (i > 0) double activity_component(int i) { assert(i > 0); return POW10(m_loggamma(i) + m_main_variables(i)); } //! Return activity of secondary species j double activity_secondary(int j) { return POW10(m_loggamma(j+m_data->nb_component))*m_concentration_secondary(j); } //! return log(IAP) for mineral m double logIAP(int m); + //! return log(IAP) for mineral "label mineral" + double logIAP(const std::string& label_mineral); //! return saturation index for mineral m double saturation_index(int m) {return logIAP(m) - m_data->logk_mineral(m);} + //! return saturation index for mineral "label mineral" + double saturation_index(const std::string& label_mineral); //! return log(IAP) for mineral m double logIAP_kinetic(int m); + //! return log(IAP) for mineral "label_mineral" + double logIAP_kinetic(const std::string& label_mineral); //! return saturation index for mineral m double saturation_index_kinetic(int m) {return logIAP(m) - m_data->logk_mineral_kinetic(m);} + //! return saturation index for mineral m + double saturation_index_kinetic(const std::string& label_mineral); private: Eigen::VectorXd m_main_variables; Eigen::VectorXd m_concentration_secondary; double m_ionic_strength; Eigen::VectorXd m_loggamma; std::shared_ptr m_data; }; } // end namespace specmicp #endif // SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP diff --git a/tests/micpsolver/test_micpsolver.hpp b/tests/micpsolver/test_micpsolver.hpp index 8595320..57ea682 100644 --- a/tests/micpsolver/test_micpsolver.hpp +++ b/tests/micpsolver/test_micpsolver.hpp @@ -1,228 +1,228 @@ /*------------------------------------------------------- - Module : tests/micpsolver - File : test_micpsolver.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "cxxtest/TestSuite.h" #include "micpsolver/micpsolver.hpp" #include "micpsolver/micpprog.hpp" #include "micpsolver/micpsolver_min.hpp" #include "utils/log.hpp" using namespace specmicp::micpsolver; class TestLinearProgram: public MiCPProg { public: //! Return the number of variables int total_variables() {return 3;} int nb_free_variables() {return 2;} int nb_complementarity_variables() {return 1;} //! Return the residual (R(X)) void get_residuals(const types::vector_t& x, types::vector_t& residual) { assert(residual.rows() == 3); residual << 1 + x(0) + x(2), 1 + 1*x(1) + x(2), -x(0) - x(1) +x(2); } //! Return the jacobian (J(x)) void get_jacobian(const types::vector_t& x, types::matrix_t& jacobian) { assert(jacobian.cols() == 3); assert(jacobian.rows() == 3); jacobian << 1, 0, 1, 0, 1, 1, -1, -1, 1; } }; class TestNonLinearProgram: public MiCPProg { public: //! Return the number of variables int total_variables() {return 3;} int nb_free_variables() {return 2;} int nb_complementarity_variables() {return 1;} //! Return the residual (R(X)) void get_residuals(const types::vector_t& x, types::vector_t& residual) { assert(residual.rows() == 3); residual << -1 + x(0)*x(0) + x(2), 1 + 1*x(1) + x(2)*x(2), -x(0) - x(1) +x(2); } //! Return the jacobian (J(x)) void get_jacobian(const types::vector_t& x, types::matrix_t& jacobian) { assert(jacobian.cols() == 3); assert(jacobian.rows() == 3); jacobian << 2*x(0), 0, 1, 0, 1, 2*x(2), -1, -1, 1; } }; /* * S. P. Dirkse and M. C. Ferris. * MCPLIB: a collection of nonlinear mixed complementarity problems. * Optimization Methods and Software, 5(4):319-345, 1995. * */ class TestKojimaProgram: public MiCPProg { public: //! Return the number of variables int total_variables() {return 4;} int nb_free_variables() {return 0;} int nb_complementarity_variables() {return 4;} //! Return the residual (R(X)) void get_residuals(const types::vector_t& x, types::vector_t& residual) { assert(residual.rows() == 4); residual << 3*x(0)*x(0)+2*x(0)*x(1)+2*x(1)*x(1)+x(2)+3*x(3)-6, 2*x(0)*x(0)+x(1)*x(1)+x(0)+10*x(2) +2*x(3)-2, 3*x(0)+x(0)*x(1)+2*x(1)*x(1)+2*x(2)+9*x(3)-9, x(0)*x(0)+3*x(1)*x(1)+2*x(2)+3*x(3)-3 ; } //! Return the jacobian (J(x)) void get_jacobian(types::vector_t& x, types::matrix_t& jacobian) { assert(jacobian.cols() == 4); assert(jacobian.rows() == 4); const int neq = total_variables(); Eigen::VectorXd res(total_variables()); Eigen::VectorXd perturbed_res(total_variables()); get_residuals(x, res); for (int j=0; j ptrprog = std::make_shared(); - MiCPSolver solver(ptrprog); + MiCPSolver solver(ptrprog); Eigen::VectorXd x(3); x << -1.5, -2 , 0; MiCPSolverReturnCode ret = solver.solve(x); TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); } void test_linear_program_10() { std::shared_ptr ptrprog = std::make_shared(); - MiCPSolver solver(ptrprog); + MiCPSolver solver(ptrprog); Eigen::VectorXd x(3); x << -1.5, -2 , 10; MiCPSolverReturnCode ret = solver.solve(x); TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); } void test_nonlinear_program_0() { std::shared_ptr ptrprog = std::make_shared(); - MiCPSolver solver(ptrprog); + MiCPSolver solver(ptrprog); Eigen::VectorXd x(3); x << -1.5, -2 , 0; MiCPSolverReturnCode ret = solver.solve(x); std::cout << x << std::endl; TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); } void test_nonlinear_program_5() { std::shared_ptr ptrprog = std::make_shared(); - MiCPSolver solver(ptrprog); + MiCPSolver solver(ptrprog); Eigen::VectorXd x(3); x << -1.5, -2 , 5; MiCPSolverReturnCode ret = solver.solve(x); std::cout << x << std::endl; TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); } void test_kojima_program() { std::shared_ptr ptrprog = std::make_shared(); - MiCPSolver solver(ptrprog); + MiCPSolver solver(ptrprog); Eigen::VectorXd x(4); x << 0.9, 0.1 , 2.9, 0.1; MiCPSolverReturnCode ret = solver.solve(x); std::cout << x << std::endl; TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); } // void test_kojima_program_min() // { // std::shared_ptr ptrprog = std::make_shared(); // MiCPSolver solver(ptrprog); // Eigen::VectorXd x(4); // x << 0.9, 0.1 , 2.9, 0.1; // MiCPSolverReturnCode ret = solver.solve(x); // std::cout << x << std::endl; // TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized); // //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8); // //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8); // //TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8); // } }; diff --git a/tests/specmicp/rate_nicoleau.cpp b/tests/specmicp/rate_nicoleau.cpp index b0c0fca..5ab38b7 100644 --- a/tests/specmicp/rate_nicoleau.cpp +++ b/tests/specmicp/rate_nicoleau.cpp @@ -1,116 +1,116 @@ #include #include #include #include #include #include // Compute the supersaturation of C3S to obtain the kinetic rate curve r=r(SI) // Data from Nicoleau et al. (2013) // // Nicoleau, L., Nonat, A., Perrey, D.: // The di- and tricalcium silicate dissolutions // Cement and Concrete Research 47(0), 14–30, 2013 const int nb_cells = 9; const double Mc3s = 228.3; using input_data_t = std::array; void read_csv_file(const std::string& filepath, std::vector& vector_data) { std::ifstream datafile(filepath); if (not datafile.is_open()) { throw std::invalid_argument("Input file cannot be opened : " + filepath); } std::string buffer; getline(datafile, buffer); // first line is headers while (datafile.good()) { input_data_t data; for (int i= 0; i data = database.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); database.swap_components(swapping); std::shared_ptr model = std::make_shared(); double m_water = 1.0; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, 0)}, {"Na[+]", specmicp::reaction_amount_t(m_naoh+m_nacl, 0)}, {"Cl[-]", specmicp::reaction_amount_t(2*m_cacl2, 0)}, {"Ca[2+]", specmicp::reaction_amount_t(m_cacl2, 0)}, {"HO[-]", specmicp::reaction_amount_t(m_naoh, 0)} }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0)}, {"Portlandite", specmicp::reaction_amount_t(m_caoh2, 0)}, }; specmicp::ReactionPathDriver driver(model, data); driver.dissolve_to_components(); Eigen::VectorXd x(data->nb_component+data->nb_mineral); x(0) = m_water; x.block(1, 0, data->nb_component-1, 1).setConstant(-2); x.block(data->nb_component, 0, data->nb_mineral, 1).setConstant(0.); specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) { std::cout << "Error : problem not solved : return code " << (int) perf.return_code << std::endl; } - return driver.get_logIAP_mineral("C3S", x); + return driver.get_current_solution().logIAP_kinetic("C3S"); } int main() { std::vector input; //read_csv_file("data_test/data_nicoleau_c3sm.csv", input); //read_csv_file("data_test/data_nicoleau_c3st1.csv", input); read_csv_file("data_test/data_nicoleau_c3st2.csv", input); // for (auto it=input.begin(); it!=input.end(); ++it) // { // for (auto its=it->begin(); its!=it->end(); ++its) // { // std::cout << *its << "\t"; // } // std::cout << std::endl; // } for (auto it=input.begin(); it!=input.end(); ++it) { std::cout << get_sursaturation(*it) << "\t" << (*it)[7] << "\t" << (*it)[8] << std::endl; } return EXIT_SUCCESS; }