diff --git a/src/specmicp/equilibrium_data.cpp b/src/specmicp/equilibrium_data.cpp index e0b8d7a..923948b 100644 --- a/src/specmicp/equilibrium_data.cpp +++ b/src/specmicp/equilibrium_data.cpp @@ -1,174 +1,173 @@ /*------------------------------------------------------- - Module : specmicp - File : equilibrium_data.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 "equilibrium_data.hpp" #include "database/module.hpp" #include "physics/constants.hpp" #define POW10(m) std::pow(10.0, m) namespace specmicp { scalar_t EquilibriumState::logIAP(index_t m) const { scalar_t logiap = 0.0; for (index_t i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral(m, i); return logiap; } scalar_t EquilibriumState::logIAP(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); index_t idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP(idm); } - scalar_t EquilibriumState::saturation_index(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); index_t idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP(idm) - m_data->logk_mineral(idm); } scalar_t EquilibriumState::logIAP_kinetic(index_t m) const { scalar_t logiap = 0.0; for (index_t i=1; inb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral_kinetic(m, i); return logiap; } scalar_t EquilibriumState::logIAP_kinetic(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); int idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals_kinetic); return logIAP_kinetic(idm); } scalar_t EquilibriumState::saturation_index_kinetic(const std::string& label_mineral) const { database::DatabaseModule namer(m_data); index_t idm = namer.safe_label_to_id(label_mineral, m_data->labels_minerals); return logIAP_kinetic(idm) - m_data->logk_mineral_kinetic(idm); } scalar_t EquilibriumState::total_aqueous_concentration_component(index_t i) const { scalar_t totconc; if (i > 0) totconc = molality_component(i); // aqueous component else totconc = 1.0/m_data->molar_mass_basis_si(0); // water totconc += m_data->nu_aqueous.col(i).dot(m_concentration_secondary); return totconc; } void EquilibriumState::total_aqueous_concentrations(Eigen::VectorXd& totconc) const { assert(totconc.rows() == m_data->nb_component); for (index_t i: m_data->range_component()) totconc(i) = total_aqueous_concentration_component(i); } void EquilibriumState::total_concentrations(Vector& total_conc) const { total_conc.resize(m_data->nb_component); total_conc.setZero(); // Aqueous concentration : m_w * B_i total_aqueous_concentrations(total_conc); total_conc = mass_water()*total_conc; // Minerals : nu_{mi} n_m for (index_t mineral: m_data->range_mineral()) { if (moles_mineral(mineral) == 0.0) continue; for (index_t component: m_data->range_component()) { total_conc(component) += m_data->nu_mineral(mineral, component)*moles_mineral(mineral); } } } scalar_t EquilibriumState::pH() const { // search where the information about "H[+]" is stored // (i.e. was the component H[+] or HO[-] ?) database::DatabaseModule namer(m_data); scalar_t id = namer.component_label_to_id("H[+]"); if (id != no_species) return -std::log10(activity_component(id)); else { id = namer.aqueous_label_to_id("H[+]"); assert(id != no_species); return -std::log10(activity_secondary(id)); } } scalar_t EquilibriumState::mass_mineral(index_t m) const { return m_data->molar_mass_mineral(m)*m_main_variables(m_data->nb_component+m); } scalar_t EquilibriumState::condensed_phase_volume() const { scalar_t volume = mass_water()/constants::water_density_25; for (index_t m=0; mnb_mineral; ++m) { volume += moles_mineral(m)*m_data->molar_volume_mineral(m); } return volume; } scalar_t EquilibriumState::gas_phase_volume() const { assert(is_fixed_volume()); return (get_volume() - condensed_phase_volume()); } scalar_t EquilibriumState::volume_minerals() const { scalar_t volume = 0; for (index_t m=0; mnb_mineral; ++m) { volume += moles_mineral(m)*m_data->molar_volume_mineral(m); } return volume; } void EquilibriumState::scale_condensed(scalar_t scale) { m_main_variables(0) *= scale; m_main_variables.segment(m_data->nb_component, m_data->nb_mineral) *= scale; } } // end namespace specmicp diff --git a/src/specmicp/reduced_system_solver.cpp b/src/specmicp/reduced_system_solver.cpp index c5cc59a..94edb84 100644 --- a/src/specmicp/reduced_system_solver.cpp +++ b/src/specmicp/reduced_system_solver.cpp @@ -1,279 +1,278 @@ /*------------------------------------------------------- - Module : specmicp - File : reduced_system_solver - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 #include "reduced_system_solver.hpp" #include "micpsolver/micpsolver.hpp" #include "micpsolver/micpsolverold.hpp" //#include "micpsolver/micpsolver_min.hpp" #include "equilibrium_data.hpp" #include "simulation_box.hpp" -#include - namespace specmicp { ReducedSystemSolver::ReducedSystemSolver( RawDatabasePtr data, const Vector& totconc, const SimulationBox& simulbox, ReducedSystemSolverOptions options ): m_options(options), m_data(data), m_system(std::make_shared( data, totconc, simulbox, options.conservation_water, options.charge_keeper )), m_var(Vector::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver( RawDatabasePtr data, const Vector& totconc, const SimulationBox& simulbox, const EquilibriumState& previous_solution, ReducedSystemSolverOptions options ): m_options(options), m_data(data), m_system(std::make_shared( data, totconc, simulbox, previous_solution, options.conservation_water, options.charge_keeper )), m_var(Vector::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver( RawDatabasePtr data, const Vector& totconc, ReducedSystemSolverOptions options ): m_options(options), m_data(data), m_system(std::make_shared( data, totconc, options.conservation_water, options.charge_keeper )), m_var(Vector::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver( RawDatabasePtr data, const Vector& totconc, const EquilibriumState& previous_solution, ReducedSystemSolverOptions options ): m_options(options), m_data(data), m_system(std::make_shared( data, totconc, SimulationBox(), previous_solution, options.conservation_water, options.charge_keeper )), m_var(Vector::Zero(data->nb_component+data->nb_mineral)) {} EquilibriumState ReducedSystemSolver::get_solution(const Eigen::VectorXd& x) { set_true_variable_vector(x); return m_system->get_solution(x, m_var); } micpsolver::MiCPPerformance ReducedSystemSolver::solve(Vector& x, bool init) { if (init) { if (get_options().ncp_function == micpsolver::NCPfunction::min) m_system->reasonable_starting_guess(x, true); else m_system->reasonable_starting_guess(x, false); } set_true_variable_vector(x); micpsolver::MiCPPerformance perf = solve(); if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized and get_options().allow_restart) { WARNING << "Failed to solve the system ! - We shake it up and start again"; const micpsolver::NCPfunction save_ncp = get_options().ncp_function; - const int save_penalization_factor = get_options().solver_options.penalization_factor; + const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor; if (save_ncp == micpsolver::NCPfunction::min) { get_options().ncp_function = micpsolver::NCPfunction::penalizedFB; m_system->reasonable_restarting_guess(x); } else { if (save_penalization_factor == 1) get_options().solver_options.penalization_factor = 0.8; set_return_vector(x); m_system->reasonable_restarting_guess(x); } set_true_variable_vector(x); micpsolver::MiCPPerformance perf2 = solve(); get_options().ncp_function = save_ncp; get_options().solver_options.penalization_factor = save_penalization_factor; perf += perf2; } if (perf.return_code > micpsolver::MiCPSolverReturnCode::NotConvergedYet) set_return_vector(x); return perf; } void ReducedSystemSolver::set_true_variable_vector(const Vector& x) { const std::vector& non_active_component = m_system->get_non_active_component(); if (non_active_component.size() == 0) { // we still copy the data, if we failed to solve the problem, we can restart if (m_options.conservation_water) m_var = x; // direct copy else m_var = x.block(1, 0, x.rows()-1, 1); for (int i=0; irange_aqueous_component()) { auto it = std::find(non_active_component.begin(), non_active_component.end(),i); if (it != non_active_component.end()) continue; m_var(new_i) = x(i); ++new_i; } for (index_t m: m_data->range_mineral()) { for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it) { if (m_data->nu_mineral(m, *it) != 0) continue; m_var(new_i) = x(m_data->nb_component+m); ++new_i; } } m_var.conservativeResize(new_i); } m_system->reduce_mineral_problem(m_var); } micpsolver::MiCPPerformance ReducedSystemSolver::solve() { if (get_options().ncp_function == micpsolver::NCPfunction::penalizedFB) { // the default micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_performance(); } else { micpsolver::MiCPSolverOLD solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_performance(); } } void ReducedSystemSolver::set_return_vector(Vector& x) { const std::vector& non_active_component = m_system->get_non_active_component(); if (non_active_component.size() == 0) { if (m_options.conservation_water) x = m_var; //direct copy else x.block(1, 0, x.rows()-1, 1) = m_var; // at that point we should have the correct solution } else { uindex_t new_i = 0; if (m_options.conservation_water) { x(0) = m_var(new_i); ++new_i; } for (index_t i: m_data->range_aqueous_component()) { auto it = std::find(non_active_component.begin(), non_active_component.end(),i); if (it != non_active_component.end()) { x(i) = -HUGE_VAL; continue; } x(i) = m_var(new_i) ; ++new_i; } for (index_t m: m_data->range_mineral()) { for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it) { if (m_data->nu_mineral(m, *it) != 0) { x(m_data->nb_component+m) = 0; continue; } x(m_data->nb_component+m) =m_var(new_i); ++new_i; } } } m_system->reset_mineral_system(m_var, x); } } // end namespace specmicp diff --git a/tests/specmicp/carboalu.cpp b/tests/specmicp/carboalu.cpp index c1088d8..c55fca4 100644 --- a/tests/specmicp/carboalu.cpp +++ b/tests/specmicp/carboalu.cpp @@ -1,205 +1,206 @@ /*------------------------------------------------------- - Module : tests/ - File : carboalu.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 -#include "specmicp/reaction_path.hpp" -#include "utils/log.hpp" #include +#include "utils/log.hpp" +#include "specmicp/reaction_path.hpp" + double mult = 1.0; double m_c3s = mult*0.6; double m_c2s = mult*0.2; double m_c3a = mult*0.10; double m_gypsum = mult*0.10; double wc = 0.8; class Timer { public: Timer() : beg_(clock_::now()) {} void reset() { beg_ = clock_::now(); } double elapsed() const { return std::chrono::duration_cast (clock_::now() - beg_).count(); } private: typedef std::chrono::high_resolution_clock clock_; typedef std::chrono::duration > second_; std::chrono::time_point beg_; }; struct Params { specmicp::micpsolver::NCPfunction ncpfunction; double penalizationfactor; bool scaling; }; const Params CCKParam = {specmicp::micpsolver::NCPfunction::penalizedFB, 0.8, false}; const Params FBParam = {specmicp::micpsolver::NCPfunction::penalizedFB, 1.0, false}; const Params minParam = {specmicp::micpsolver::NCPfunction::min, 1.0, true}; void solve_carboalu(double delta, Params param, bool output=true, bool coldstart=false) { specmicp::database::Database database("../data/cemdata_specmicp.js"); std::shared_ptr data = database.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Al[3+]","Al(OH)4[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); database.swap_components(swapping); std::shared_ptr model = std::make_shared(); double delta_h2co3 =delta; double delta_sio2 = 0.00; double m_water = wc*1e-3*(m_c3s*(3*56.08+60.08)+m_c2s*(2*56.06+60.08)+m_c3a*(3*56.08+101.96)+m_gypsum*(56.08+80.06+2*18.02)); //std::cout << m_water << std::endl; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/data->molar_mass_basis_si(0), delta_h2co3+delta_sio2)}, {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)}, {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3-delta_sio2)}, {"SiO(OH)3[-]", specmicp::reaction_amount_t(0, delta_sio2)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)} }; - model->database_path = "data/cemdata_specmicp.js"; + model->database_path = "../data/cemdata_specmicp.js"; model->minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Calcite", "Al(OH)3,am", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", "Straetlingite", "Gypsum", "Ettringite", "Thaumasite" }; Eigen::VectorXd x; specmicp::ReactionPathDriver driver(model, data, x); Eigen::VectorXd totaq(data->nb_component); double totiter = 0; double totfact = 0; driver.get_options().solver_options.penalization_factor = param.penalizationfactor; driver.get_options().ncp_function = param.ncpfunction; driver.get_options().solver_options.use_scaling = param.scaling; driver.get_options().solver_options.max_factorization_step = 4; driver.get_options().solver_options.factor_gradient_search_direction = 200; driver.get_options().solver_options.max_iter = 50; driver.get_options().solver_options.maxstep = 100; + driver.get_options().solver_options.factor_descent_condition = 1e-6; driver.get_options().allow_restart = true; driver.get_options().charge_keeper = database.component_label_to_id("HO[-]"); if (output) { - std::cout << "Reaction_path return_code nb_iter nb_fact pH"; - for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it) - { + std::cout << "Reaction_path return_code nb_iter nb_fact pH"; + for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it) + { std::cout << " " << *it; - } - for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it) - { + } + for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it) + { std::cout << " " << *it; - } - std::cout << std::endl; + } + std::cout << std::endl; } const int nb_step = 1 + 2.5/delta; for (int i=0; inb_component-1,1).transpose(); - for (int m=0; m < data->nb_mineral; ++m) std::cout << " " << solution.moles_mineral(m); - std::cout << std::endl; + std::cout << i*(delta_sio2+delta_h2co3) << " " << (int) perf.return_code << " " + << perf.nb_iterations << " " << perf.nb_factorization << " " + << solution.pH() << " " + << solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose(); + for (int m=0; m < data->nb_mineral; ++m) std::cout << " " << solution.moles_mineral(m); + std::cout << std::endl; } totiter += perf.nb_iterations; totfact += perf.nb_factorization; //data->stability_mineral(database.mineral_label_to_id("SiO2,am")) = specmicp::database::MineralStabilityClass::cannot_dissolve; } if (output) { std::cout << "Average iterations : " << totiter/nb_step << std::endl; std::cout << "Average factorization : " << totfact/nb_step << std::endl; } } int main() { - specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::logger::ErrFile::stream() = &std::cerr; int nbloop = 7500; double delta = 0.01; Timer tmr; solve_carboalu(delta, CCKParam, true); //for (int i=0; i