diff --git a/src/specmicp/problem_solver/step_solver.cpp b/src/specmicp/problem_solver/step_solver.cpp index 22712f7..1446dc4 100644 --- a/src/specmicp/problem_solver/step_solver.cpp +++ b/src/specmicp/problem_solver/step_solver.cpp @@ -1,313 +1,326 @@ /* ============================================================================= 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 "step_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/problem_solver/smart_solver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/problem_solver/reactant_box.hpp" #include "specmicp_common/log.hpp" #include #include namespace specmicp { // =========== // // StepProblem // // =========== // struct StepProblem::StepProblemImpl { StepProblemImpl() {} stepper_next_step_f m_constraints_updater; stepper_stop_condition_f m_stopper; RawDatabasePtr m_data; units::UnitsSet m_units; AdimensionalSystemConstraints m_constraints; AdimensionalSystemSolutionExtractor get_extractor( const AdimensionalSystemSolution * const solution ) { return AdimensionalSystemSolutionExtractor(*solution, m_data, m_units); } }; StepProblem::StepProblem(): m_impl(utils::make_pimpl()) {} StepProblem::StepProblem( stepper_next_step_f updater, stepper_stop_condition_f stopper, AdimensionalSystemConstraints&& constraints ): m_impl(utils::make_pimpl()) { m_impl->m_constraints_updater = updater; m_impl->m_stopper = stopper; m_impl->m_constraints = std::move(constraints); } StepProblem::StepProblem(StepProblem &&other): m_impl(std::move(other.m_impl)) { } StepProblem::~StepProblem() = default; void StepProblem::set_constraint_updater(stepper_next_step_f updater) { m_impl->m_constraints_updater = updater; } void StepProblem::set_stop_condition(stepper_stop_condition_f stopper) { m_impl->m_stopper = stopper; } void StepProblem::set_stop_condition_step(uindex_t stop_step) { set_stop_condition([stop_step]( specmicp::uindex_t cnt, const specmicp::AdimensionalSystemSolutionExtractor& _) { auto status = specmicp::StepProblemStatus::OK; if (cnt >= stop_step) {status = specmicp::StepProblemStatus::Stop;} return status; }); } void StepProblem::set_initial_constraints( ReactantBox& reactant_box, bool modify_db ) { m_impl->m_constraints = reactant_box.get_constraints(modify_db); m_impl->m_data = reactant_box.get_database(); m_impl->m_units = reactant_box.get_units(); } bool StepProblem::is_valid() { return ( m_impl->m_constraints_updater != nullptr and m_impl->m_stopper != nullptr); } AdimensionalSystemConstraints StepProblem::get_constraints() const { return m_impl->m_constraints; } StepProblemStatus StepProblem::update_constraints( uindex_t cnt, const AdimensionalSystemSolution * const solution ) { if (m_impl->m_constraints_updater == nullptr) { ERROR << "Bad initialization : no constraints updater in step problem."; } auto status = m_impl->m_constraints_updater( cnt, m_impl->get_extractor(solution), m_impl->m_constraints ); return status; } StepProblemStatus StepProblem::check_stop_condition( uindex_t cnt, const AdimensionalSystemSolution * const solution) { if (m_impl->m_stopper == nullptr) { ERROR << "Bad initialization : no stopper in step problem."; } return m_impl->m_stopper(cnt, m_impl->get_extractor(solution) ); } std::shared_ptr StepProblem::get_database() { return m_impl->m_data; } const units::UnitsSet& StepProblem::get_units() { return m_impl->m_units; } // ================= // // StepProblemRunner // // ================= // struct StepProblemRunner::StepProblemRunnerImpl { AdimensionalSystemSolverOptions m_opts; AdimensionalSystemSolution m_solution; stepper_output_f m_output {nullptr}; }; StepProblemRunner::StepProblemRunner(): m_impl(utils::make_pimpl()) { } StepProblemRunner::StepProblemRunner(StepProblemRunner&& other): m_impl(std::move(other.m_impl)) { } StepProblemRunner::~StepProblemRunner() = default; void StepProblemRunner::set_warmstart_solution(const AdimensionalSystemSolution& solution) { m_impl->m_solution = solution; } void StepProblemRunner::set_warmstart_solution(AdimensionalSystemSolution&& solution) { m_impl->m_solution = std::move(solution); } void StepProblemRunner::set_solver_options(AdimensionalSystemSolverOptions& opts) { m_impl->m_opts = opts; } AdimensionalSystemSolverOptions& StepProblemRunner::get_solver_options() { return m_impl->m_opts; } const AdimensionalSystemSolution& StepProblemRunner::get_solution() const { return m_impl->m_solution; } const AdimensionalSystemSolution * const StepProblemRunner::get_ptr_solution() const { return &m_impl->m_solution; } void StepProblemRunner::set_step_output(stepper_output_f output) { m_impl->m_output = output; } StepProblemStatus StepProblemRunner::run(StepProblem& problem) { uindex_t cnt = 0; StepProblemStatus status; while (cnt < std::numeric_limits::max()) { SmartAdimSolver solver(problem.get_database(), problem.get_constraints(), get_solver_options() ); if (get_solution().is_valid) { solver.set_warmstart(get_solution()); } bool ret = solver.solve(); if (not ret) { status = StepProblemStatus::Error; ERROR << "Failed to solve step " << cnt << " in step problem"; break; } // save solution m_impl->m_solution = solver.get_solution(); // output if asked if (m_impl->m_output != nullptr) { const AdimensionalSystemSolution& solref = m_impl->m_solution; status = m_impl->m_output( cnt, AdimensionalSystemSolutionExtractor( solref, problem.get_database(), problem.get_units()) ); if (status == StepProblemStatus::Error) { ERROR << "User detected problem in step " << cnt << " while updating contraints "; break; } } ++cnt; status = problem.check_stop_condition(cnt, get_ptr_solution()); if (status == StepProblemStatus::Stop) break; else if (status == StepProblemStatus::Error) { ERROR << "User detected problem in step " << cnt; break; } status = problem.update_constraints(cnt, get_ptr_solution()); if (status == StepProblemStatus::Error) { ERROR << "User detected problem in step " << cnt << "while updating contraints "; break; } } if (status == StepProblemStatus::Stop) status = StepProblemStatus::OK; else status = StepProblemStatus::Error; return status; } +// Functors +// ========= + + +StepProblemStatus StepTotalConcentration::operator() ( + uindex_t cnt, + const AdimensionalSystemSolutionExtractor& _, + AdimensionalSystemConstraints& constraints) +{ + constraints.total_concentrations += m_update; + return StepProblemStatus::OK; +} + } // end namespace specmicp diff --git a/src/specmicp/problem_solver/step_solver.hpp b/src/specmicp/problem_solver/step_solver.hpp index ed50a5e..771dc84 100644 --- a/src/specmicp/problem_solver/step_solver.hpp +++ b/src/specmicp/problem_solver/step_solver.hpp @@ -1,211 +1,248 @@ /* ============================================================================= 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_STEPSOLVER_HPP #define SPECMICP_SPECMICP_STEPSOLVER_HPP /*! * \file specmicp/problem_solver/step_solver.hpp * \brief Solve a series of problem + * + * It provides a common interface to solve a succession of similar problem. */ #include "specmicp_common/types.hpp" #include "specmicp_common/pimpl_ptr.hpp" #include #include namespace specmicp { struct AdimensionalSystemConstraints; struct AdimensionalSystemSolution; class AdimensionalSystemSolutionExtractor; struct AdimensionalSystemSolverOptions; class ReactantBox; namespace database { struct DataContainer; } // end namespace database namespace units { struct UnitsSet; } // end namespace units //! \brief The current status of the computation enum class StepProblemStatus { OK, //!< Continue the computation Stop, //!< The computation finished correctly Error //!< An error occured, the computation was aborted }; //! \brief A function of this type check when the computation must be stopped using stepper_stop_condition_f = std::function; //! \brief A function of this type computes the next set of constraints //! //! return : status of the computation //! param : //! - uindex_t [in] : the current step id //! - const AdimensionalSystemSolutionExtractor& [in] : the previous solution //! - AdimensionalSystemConstraints& [out] : the set of constraints to update using stepper_next_step_f = std::function; //! \brief A function of this type output the solution at each step using stepper_output_f = std::function; //! \brief A step problem //! //! A step problem is a series of computations, each only slightly different. class SPECMICP_DLL_PUBLIC StepProblem { public: StepProblem(); StepProblem(stepper_next_step_f updater, stepper_stop_condition_f stopper, AdimensionalSystemConstraints&& constraints ); StepProblem(StepProblem&& other); ~StepProblem(); // Setter // ====== //! \brief Set how to compute the next set of constraints void set_constraint_updater(stepper_next_step_f updater); //! \brief Set how to end the computation void set_stop_condition(stepper_stop_condition_f stopper); //! \brief Stop the problem when cnt == max_step void set_stop_condition_step(uindex_t stop_step); //! \brief Set initial constraints void set_initial_constraints(ReactantBox& reactant_box, bool modify_db=true); //! \brief Return true if the problem is correctly initialized bool is_valid(); // Getter // ====== AdimensionalSystemConstraints get_constraints() const; //! \brief Update constraints StepProblemStatus update_constraints( uindex_t cnt, const AdimensionalSystemSolution * const solution ); //! \brief Check if the computation is finished //! //! \param cnt current step //! \param solution current solution StepProblemStatus check_stop_condition( uindex_t cnt, const AdimensionalSystemSolution * const solution ); //! \brief Return the database std::shared_ptr get_database(); //! \brief Return the units const units::UnitsSet& get_units(); private: StepProblem(const StepProblem& other) = delete; StepProblem& operator=(const StepProblem& other) = delete; struct SPECMICP_DLL_LOCAL StepProblemImpl; utils::pimpl_ptr m_impl; }; + +//! \brief Run a step problem +//! +//! This is the 'loop' algorithm which solves every step. class SPECMICP_DLL_PUBLIC StepProblemRunner { public: StepProblemRunner(); StepProblemRunner(StepProblemRunner&& other); ~StepProblemRunner(); //! \brief Set an initial solution to warmastart the computation (optional) void set_warmstart_solution(const AdimensionalSystemSolution& solution); void set_warmstart_solution(AdimensionalSystemSolution&& solution); //! \brief Set the specmicp solver options void set_solver_options(AdimensionalSystemSolverOptions& opts); //! \brief Set the callback to save the solution at each step void set_step_output(stepper_output_f output); //! \brief Return a reference to the solver options AdimensionalSystemSolverOptions& get_solver_options(); //! \brief Return the solution const AdimensionalSystemSolution& get_solution() const; //! \brief Solve the problem StepProblemStatus run(StepProblem& problem); private: StepProblemRunner(const StepProblem& other) = delete; StepProblemRunner& operator=(const StepProblem& other) = delete; struct SPECMICP_DLL_LOCAL StepProblemRunnerImpl; utils::pimpl_ptr m_impl; //! \brief Return the solution const AdimensionalSystemSolution * const get_ptr_solution() const; }; + +// Simple functors to use as constraint updaters + +//! \brief Update the total concentration at each step +//! +//! \param update_concentration update to the total concentrations +//! +//! The update must be a vector Nc x 1 +class StepTotalConcentration +{ +public: + StepTotalConcentration( + const Eigen::Ref& update_concentration + ): + m_update(update_concentration) + {} + StepTotalConcentration( + Vector&& update_concentration + ): + m_update(update_concentration) + {} + + StepProblemStatus operator() ( + uindex_t cnt, + const AdimensionalSystemSolutionExtractor& _, + AdimensionalSystemConstraints& constraints); + +private: + Vector m_update; +}; + } // end namespace specmicp #endif // SPECMICP_SPECMICP_STEPSOLVER_HPP diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index bc65550..4d9b68f 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,530 +1,534 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/problem_solver/reactant_box.hpp" #include "specmicp/problem_solver/smart_solver.hpp" #include "specmicp/problem_solver/step_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/physics/laws.hpp" #include TEST_CASE("Reactant Box", "[units],[init],[formulation]") { SECTION("Units setting - SI") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); auto raw_data= thedatabase.get_database(); auto units_set = specmicp::units::SI_units; specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_aqueous_species("CaCl2", 5.0, "kg"); reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg"); auto vector = reactant_box.get_total_concentration(false); auto mass_sol = 0.5 * specmicp::laws::density_water(units_set); CHECK(vector(0) == Approx( mass_sol / raw_data->molar_mass_basis(0, units_set) )); auto id_ca = raw_data->get_id_component("Ca[2+]"); auto id_na = raw_data->get_id_component("Na[+]"); auto id_cl = raw_data->get_id_component("Cl[-]"); auto id_cacl2 = raw_data->get_id_compound("CaCl2"); CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2))); CHECK(vector(id_na) == Approx(0.5*mass_sol)); CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca))); } SECTION("Units setting - mol L") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); auto raw_data= thedatabase.get_database(); auto units_set = specmicp::units::SI_units; units_set.length = specmicp::units::LengthUnit::decimeter; specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(0.5, "L"); reactant_box.add_aqueous_species("CaCl2", 5.0, "kg"); reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg"); auto vector = reactant_box.get_total_concentration(false); auto mass_sol = 0.5 * specmicp::laws::density_water(units_set); CHECK(vector(0) == Approx( mass_sol / raw_data->molar_mass_basis(0, units_set) )); auto id_ca = raw_data->get_id_component("Ca[2+]"); auto id_na = raw_data->get_id_component("Na[+]"); auto id_cl = raw_data->get_id_component("Cl[-]"); auto id_cacl2 = raw_data->get_id_compound("CaCl2"); CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2))); CHECK(vector(id_na) == Approx(0.5*mass_sol)); CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca))); } } TEST_CASE("Solving adimensional system with problem setting", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Solving problem with dissolver - simpler problem") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; formulation.mass_solution = 0.8; formulation.amount_minerals = {{"Portlandite", 10.0}}; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); total_concentrations *= 1e3; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); solver.get_options().solver_options.maxstep = 50.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.factor_descent_condition = 1e-10; solver.get_options().solver_options.factor_gradient_search_direction = 200; specmicp::Vector x; solver.initialize_variables(x, 0.8, -2); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.802367).epsilon(1e-4)); CHECK(extr.volume_fraction_mineral(raw_data->get_id_mineral("Portlandite")) == Approx(0.32947).epsilon(1e-2)); } SECTION("Solving problem with dissolver") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; formulation.mass_solution = 0.8; formulation.amount_minerals = {{"C3S", 0.7}, {"C2S", 0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); total_concentrations *= 1e3; total_concentrations(2) = 500; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); solver.get_options().solver_options.maxstep = 20.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.factor_descent_condition = 1e-10; solver.get_options().solver_options.factor_gradient_search_direction = 200; specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral()); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); } SECTION("Solving problem with reactant box") { std::cout << "\n Solving with reactant box \n ------------------------ \n"; std::cout << "== with SI units == \n"; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"HCO3[-]", "CO2"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); reactant_box.set_charge_keeper("HO[-]"); specmicp::AdimensionalSystemConstraints constraints = reactant_box.get_constraints(true); specmicp::Vector total_concentrations = constraints.total_concentrations; specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.5, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(solution,raw_data, specmicp::units::SI_units); std::cout << "== with non SI units : mmol/dm == \n"; specmicp::units::UnitsSet dmmol; dmmol.length = specmicp::units::LengthUnit::decimeter; dmmol.quantity = specmicp::units::QuantityUnit::millimoles; reactant_box.set_units(dmmol); specmicp::AdimensionalSystemConstraints constraints2 = reactant_box.get_constraints(false); specmicp::Vector total_concentrations2 = constraints2.total_concentrations; specmicp::AdimensionalSystemSolver solver2(raw_data, constraints2); solver2.get_options().units_set = dmmol; solver2.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::Vector x2; solver2.initialize_variables(x2, 0.5, -6); perf = solver2.solve(x2); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; specmicp::AdimensionalSystemSolution solution2 = solver2.get_raw_solution(x); CHECK(solution.main_variables(0) == Approx(solution2.main_variables(0)).epsilon(1e-6)); CHECK(solution.main_variables(2) == Approx(solution2.main_variables(2)).epsilon(1e-6)); CHECK(solution.main_variables(5) == Approx(solution2.main_variables(5)).epsilon(1e-6)); CHECK(solution.ionic_strength == Approx(solution2.ionic_strength).epsilon(1e-6)); specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, dmmol); for (auto idc: raw_data->range_component()) { if (idc == specmicp::database::electron_id) continue; // std::cout << "Id " << idc << " - " << raw_data->get_label_component(idc) << std::endl; // std::cout << total_concentrations(idc) << " - " // << total_concentrations2(idc) << " - " // << extr.total_concentration(idc) << " ~ " // << extr.volume_fraction_gas_phase() * extr.total_gaseous_concentration(idc) << " # " // << extr.fugacity_gas(0) << " - " // << extr2.total_concentration(idc) // << std::endl; CHECK(extr.total_concentration(idc) == Approx(total_concentrations(idc)).epsilon(1e-8)); CHECK(extr2.total_concentration(idc) == Approx(total_concentrations2(idc)).epsilon(1e-8)); CHECK(extr.total_concentration(idc) == Approx(extr2.total_concentration(idc)).epsilon(1e-8)); } std::cout << "== with non SI units : mmol/cm == \n"; specmicp::units::UnitsSet cmmol; cmmol.length = specmicp::units::LengthUnit::centimeter; cmmol.quantity = specmicp::units::QuantityUnit::millimoles; reactant_box.set_units(cmmol); specmicp::AdimensionalSystemConstraints constraints3 = reactant_box.get_constraints(false); specmicp::Vector total_concentrations3 = constraints3.total_concentrations; specmicp::AdimensionalSystemSolver solver3(raw_data, constraints3); solver3.get_options().units_set = cmmol; solver3.get_options().solver_options.set_tolerance(1e-8, 1e-12); specmicp::Vector x3; solver3.initialize_variables(x3, 0.5, -6); perf = solver3.solve(x3); for (auto idc: raw_data->range_component()) { if (idc == specmicp::database::electron_id) continue; CHECK(total_concentrations(idc) == Approx(1e3*total_concentrations3(idc)).epsilon(1e-8)); } REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; specmicp::AdimensionalSystemSolution solution3 = solver3.get_raw_solution(x); CHECK(solution.main_variables(0) == Approx(solution3.main_variables(0)).epsilon(1e-6)); CHECK(solution.main_variables(2) == Approx(solution3.main_variables(2)).epsilon(1e-6)); CHECK(solution.main_variables(5) == Approx(solution3.main_variables(5)).epsilon(1e-6)); CHECK(solution.ionic_strength == Approx(solution3.ionic_strength).epsilon(1e-6)); specmicp::AdimensionalSystemSolutionExtractor extr3(solution3, raw_data, cmmol); for (auto idc: raw_data->range_component()) { CHECK(extr3.total_concentration(idc) == Approx(total_concentrations3(idc)).epsilon(1e-8)); } } SECTION("Solving problem with reactant box and smart solver") { std::cout << "\n Solving with reactant box and smart solver \n ------------------------ \n"; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"HCO3[-]", "CO2"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); reactant_box.set_charge_keeper("HO[-]"); specmicp::SmartAdimSolver solver(raw_data, reactant_box); solver.set_init_volfrac_water(0.5); solver.set_init_molality(1e-6); solver.solve(); REQUIRE(solver.is_solved() == true); } SECTION("Reactant box fixed activity") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(0.5, "L/L"); reactant_box.add_solid_phase("C3S", 0.350, "L/L"); reactant_box.add_solid_phase("C2S", 0.05, "L/L"); reactant_box.set_charge_keeper("HCO3[-]"); reactant_box.add_fixed_activity_component("H[+]", 1e-9); auto constraints = reactant_box.get_constraints(true); CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]")); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.6, -4); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.pH() == Approx(9)); } SECTION("Reactant box fixed molality") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"}, // {"H[+]", "HO[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(0.5, "L/L"); reactant_box.add_solid_phase("C3S", 0.35, "L/L"); reactant_box.add_solid_phase("C2S", 0.05, "L/L"); reactant_box.set_charge_keeper("HCO3[-]"); reactant_box.add_fixed_molality_component("H[+]", 1e-5); auto constraints = reactant_box.get_constraints(true); CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]")); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.6, -4); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.molality_component(raw_data->get_id_component("H[+]")) == Approx(1e-5)); } SECTION("Reactant box fixed fugacity") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_fixed_fugacity_gas("CO2(g)", "HCO3[-]", 1e-3); reactant_box.set_charge_keeper("H[+]"); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.8, -6); solver.get_options().solver_options.set_tolerance(1e-6, 1e-8); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.fugacity_gas(raw_data->get_id_gas("CO2(g)")) == Approx(1e-3)); } SECTION("Solving problem with stepper") { std::cout << "\n Solving with stepper \n ------------------------ \n"; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"HCO3[-]", "CO2"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); + // initial constraints specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); - //reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); reactant_box.set_charge_keeper("HO[-]"); reactant_box.add_component("CO2"); - specmicp::uindex_t max_step = 100; - std::vector pHs; - pHs.reserve(max_step); + specmicp::uindex_t max_step = 20; specmicp::StepProblem problem; problem.set_initial_constraints(reactant_box, true); specmicp::index_t id_co2 = raw_data->get_id_component("CO2"); specmicp::index_t id_ca = raw_data->get_id_component("Ca[2+]"); + // The updates to apply at each step + specmicp::Vector update; + update.setZero(raw_data->nb_component()); + update(id_co2) = (0.95*problem.get_constraints().total_concentrations(id_ca))/max_step; + + problem.set_constraint_updater(specmicp::StepTotalConcentration(update)); problem.set_stop_condition_step(max_step); - problem.set_constraint_updater( - [max_step, raw_data, id_co2, id_ca]( - specmicp::uindex_t cnt, - const specmicp::AdimensionalSystemSolutionExtractor& _, - specmicp::AdimensionalSystemConstraints& constr){ - constr.total_concentrations(id_co2) = double(cnt)/max_step*(0.95*constr.total_concentrations(id_ca)); - return specmicp::StepProblemStatus::OK; - }); + specmicp::StepProblemRunner runner; + + // output + std::vector pHs; + pHs.reserve(max_step); runner.set_step_output( [&pHs](specmicp::uindex_t cnt, const specmicp::AdimensionalSystemSolutionExtractor& extr) { pHs.push_back(extr.pH()); return specmicp::StepProblemStatus::OK; }); + + + auto status = runner.run(problem); REQUIRE(status == specmicp::StepProblemStatus::OK); CHECK(pHs.size() == max_step); for (auto it: pHs) { std::cout << it << std::endl; } - CHECK(pHs[pHs.size()-1] == Approx(9.84065).epsilon(1e-4)) + CHECK(pHs[pHs.size()-1] == Approx(9.84065).epsilon(1e-4)); } }