diff --git a/CMakeLists.txt b/CMakeLists.txt index 3e3494d..9d418e7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,239 +1,241 @@ project(specmicp) cmake_minimum_required(VERSION 2.8) ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package find_package(CxxTest) # Necessary if you want the tests ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### set(PROJECT_TEST_DIR ${PROJECT_SOURCE_DIR}/tests) set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${EIGEN3_INCLUDE_DIR}) include_directories(${PROJECT_SOURCE_DIR}) ################ Modules ########################## # -------- MiCPSolver ---------------- set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(specmisolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl + ${MICPSOLVER_DIR}/micpsolverold.hpp + ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # ------- SpecMiCP --------------- set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/thermodata.hpp ${SPECMICP_DIR}/simulation_box.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ${SPECMICP_DIR}/speciation_solver.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # -------- Database ---------------- set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/data_container.cpp ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # --------- ODE int ------------------ set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ) # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp ) # ------- Physics ---------------- set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # ------- Data ------------------- set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # ------ Documentation ----------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) endif(DOXYGEN_FOUND) # --------- Test ---------------------- # CXX Test # ======== if(CXXTEST_FOUND) set(CXXTEST_USE_PYTHON TRUE) include_directories(${CXXTEST_INCLUDE_DIR}) enable_testing() # MiCP Solver # ------------ CXXTEST_ADD_TEST(test_cond_number test_cond_number.cpp ${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp) CXXTEST_ADD_TEST(test_ncp_funtion test_ncp_function.cpp ${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp) CXXTEST_ADD_TEST(test_micpsolver test_micpsolver.cpp ${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp) # Database # -------- CXXTEST_ADD_TEST(test_data_reader test_data_reader.cpp ${PROJECT_TEST_DIR}/database/test_data_reader.hpp) target_link_libraries(test_data_reader specmicp_database jsoncpp) CXXTEST_ADD_TEST(test_data_selector test_data_selector.cpp ${PROJECT_TEST_DIR}/database/test_data_selector.hpp) target_link_libraries(test_data_selector specmicp_database jsoncpp) # Odeint # ------ CXXTEST_ADD_TEST(test_butchertableau test_butchertableau.cpp ${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp) CXXTEST_ADD_TEST(test_embeddedrungekutta test_embeddedrungekutta.cpp ${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp) endif() # Test 'Perso' # ============ add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp) target_link_libraries(thermocarbo specmicp specmicp_database jsoncpp) add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp) target_link_libraries(carboalu specmicp specmicp_database jsoncpp) add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp) target_link_libraries(alkaliactivated specmicp specmicp_database jsoncpp) add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp) target_link_libraries(rate_nicoleau specmicp specmicp_database jsoncpp) add_executable(kinetics ${PROJECT_TEST_DIR}/specmicp/kinetics.cpp) target_link_libraries(kinetics specmicp specmicp_database jsoncpp) set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv ) add_custom_target(data_nicoleau SOURCES ${DATA_TEST}) file(INSTALL ${DATA_TEST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/) diff --git a/src/micpsolver/micpsolverold.hpp b/src/micpsolver/micpsolverold.hpp new file mode 100644 index 0000000..176ec07 --- /dev/null +++ b/src/micpsolver/micpsolverold.hpp @@ -0,0 +1,302 @@ +/*------------------------------------------------------- + + - Module : micpsolver + - File : micpsolver.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#ifndef SPECMIC_MICPSOLVER_MICPSOLVEROLD_HPP +#define SPECMIC_MICPSOLVER_MICPSOLVEROLD_HPP + +#include +#include +#include "micpsolver_structs.hpp" +#include "ncp_function.hpp" + +//! \file micpsolver.hpp The MiCP solver + +namespace specmicp { +namespace micpsolver { + +//! \brief Call a NCP-function +//! +//! @tparam ncp_t the NCP function to use +//! @param a first argument +//! @param b second argument +//! @param t parameter of the NCP function +//! +template +inline double ncp_function(double a, double b, double t); + +//! \brief Reformulate the jacobian using the NCP-function ncp_t +//! +//! @tparam ncp_t the NCP-function to use +//! @param[in] neq number of equation +//! @param[in] neq_free number of free variables +//! @param[in] x the variables +//! @param[in] r the residuals +//! @param[in,out] jacobian the jacobian to reformulate +//! @param[in,out] r_reformulated the reformulated residuals +//! @param[in] t parameter of the NCP function +template +inline void reformulate_jacobian_helper( + int neq, + int neq_free, + const Eigen::VectorXd& x, + const Eigen::VectorXd& r, + Eigen::MatrixXd& jacobian, + Eigen::VectorXd& r_reformulated, + double t + ); + +//! \brief reformulate the result of the Newton system +template +inline void reformulate_result( + int neq, + int neq_free, + Eigen::VectorXd& x, + const Eigen::VectorXd& orig_r, + Eigen::VectorXd& grad_phi, + Eigen::VectorXd& update + ); + +//! \brief The MiCP Solver +//! +//! Solve +//! - \f$ G(u, v) = 0\f$ +//! - \f$ 0 \leq v \perp H(u,v) \geq 0 \f$ +//! +//! \tparam Program a subclass of MiCPProg +//! +//! References : +//! - \cite Munson2001 +//! - \cite Facchinei2003 +//! +template +class MiCPSolverOLD { + +public: + //! \brief Constructor + //! + //! @param prog smart pointer toward an instance of a Program to solve + MiCPSolverOLD(std::shared_ptr prog): + m_program(prog), + m_residuals(prog->total_variables()), + m_phi_residuals(Eigen::VectorXd::Zero(prog->total_variables())), + m_jacobian(prog->total_variables(),prog ->total_variables()), + m_max_taken(false), + m_consec_max(0) + { + } + + //! \brief Return a const reference to the options used by the algorithm + const MiCPSolverOptions& get_options() const {return m_options;} + //! \brief Return a reference to the options used by the algorithm + MiCPSolverOptions& get_options() {return m_options;} + //! \brief Set the options + void set_options(const MiCPSolverOptions& options) {m_options = options;} + + //! \brief Return the number of equations + int get_neq() const {return m_program->total_variables();} + //! \brief Return the number of equations corresponding to the free variables (size of G) + int get_neq_free() const {return m_program->nb_free_variables();} + + // Merit function + // ############## + + //! \brief Reformulation for lower bounded variable + double phi_lower_bounded(const double& x, const double& r, const double& l) const { + return ncp_function(x-l, r, get_options().penalization_factor);} + //! \brief Reformulation for a free variable + double phi_free(const double& r) const { + return r; + } + + //! \brief Compute the norm of the update + template + double norm_update(const Eigen::VectorXd& update, + const Eigen::VectorXd& solution) const { + return (update.array().abs()/(solution.array().max(1)) + ).matrix().template lpNorm

(); + } + // Residuals and jacobian + // ---------------------- + + //! \brief Compute the residuals, store it in r + //! + //! @param[in] x the variables + //! @param[out] r vector to store the residuals (must be of the same size as x) + void compute_residuals(const Eigen::VectorXd& x, Eigen::VectorXd& r) + { + m_program->get_residuals(x, r); + get_perf().nb_call_residuals += 1; + } + //! \brief Compute the residuals, use internal storage + //! + //! @param[in] x the variables + void compute_residuals(const Eigen::VectorXd& x) + { + return compute_residuals(x, m_residuals); + get_perf().nb_call_jacobian += 1; + } + //! \brief Compute the jacobian + //! + //! Assumes that the residual have been computed before + //! + //! @param[in] x the variables + void compute_jacobian(Eigen::VectorXd& x) + { + m_program->get_jacobian(x, m_jacobian); + } + //! \brief Reformulation of the residuals + //! + //! Reformulate the problem - assumes that r, the residual, has been computed before + //! + //! @param[in] x the variables + //! @param[in] r the residuals + //! @param[out] r_phi a vector of size neq, which will contain the reformulated residuals + void reformulate_residuals(const Eigen::VectorXd& x, + const Eigen::VectorXd& r, + Eigen::VectorXd& r_phi); + //! \brief Reformulation of the residuals *inplace* + //! + //! @param[in] x the variables + //! @param[in,out] r the residual, will contain the reformulated residuals + void reformulate_residuals_inplace(const Eigen::VectorXd& x, + Eigen::VectorXd &r); + //! \brief Reformulation of the jacobian + //! + //! r is the original vector of residuals (not reformulated) + void reformulate_jacobian(const Eigen::VectorXd& x + ) + { + reformulate_jacobian_helper(get_neq(), get_neq_free(), + x, m_residuals, + m_jacobian, + m_phi_residuals, + get_options().penalization_factor); + } + + //! \brief Compute the factors to scale the jacobian + //! + //! @param[in] jacobian the jacobian to scale (from the reformulated problem) + //! @param[in] residual the residuals corresponding to the jacobian + //! @param[out] rscaler scaling factors of the rows + //! @param[out] cscaler scaling factors of the columns + void scaling_jacobian(const Eigen::MatrixXd& jacobian, + const Eigen::VectorXd& residual, + Eigen::VectorXd& rscaler, + Eigen::VectorXd& cscaler); + + // Algorithm + // ######### + + //! \brief Solver the program using x as initial guess + //! + //! @param[in,out] x the initial guess, as output contains the solution (from the last iteration) + MiCPSolverReturnCode solve(Eigen::VectorXd& x); + + //! \brief Setup the residuals + //! + //! @param[in] x the current solution + void setup_residuals(const Eigen::VectorXd& x) { + compute_residuals(x); + reformulate_residuals(x, m_residuals, m_phi_residuals); + } + + //! \brief Setup the jacobian + //! + //! @param[in] x the current solution + void setup_jacobian(Eigen::VectorXd& x) { + compute_jacobian(x); + reformulate_jacobian(x); + m_grad_phi = m_jacobian.transpose() * m_phi_residuals; + } + + //! \brief Check for convergence + //! + //! @param nb_iterations the number of iterations + //! @param update the update taken at the previous iteration + //! @param solution the current solution + //! @return a MiCPSolverReturnCode describing the state of the algorithm + MiCPSolverReturnCode check_convergence(int nb_iterations, + Eigen::VectorXd& update, + Eigen::VectorXd& solution); + + //! \brief Solve the Newton system + //! + //! Assume that the Newton system has been formed + //! + //! \param[out] update the update to apply to the solution + MiCPSolverReturnCode search_direction_calculation(Eigen::VectorXd& update); + + //! \brief Solve the Newton system - does not scale the jacobian + //! + //! Assume that the Newton system has been formed + //! + //! \param[out] update the update to apply to the solution + MiCPSolverReturnCode search_direction_calculation_no_scaling(Eigen::VectorXd& update); + //! \brief Linesearch + //! + //! Nonmonotone linesearch + //! + //! References : + //! - \cite Dennis1983 + //! - \cite Munson2001 + //! + int linesearch(Eigen::VectorXd& update, Eigen::VectorXd& x); + + //! \brief Crashing + //! + //! This function improves, if possible, the initial guess + //! Reference : + //! - \cite Munson2001 + void crashing(Eigen::VectorXd& x); + + //! \brief Project variables on the feasible set + void projection(Eigen::VectorXd& x); + + //! \brief Return a const reference to an instance of MiCPPerformance + const MiCPPerformance& get_performance() {return m_perf;} + +protected: + //! \brief Return a reference to an instance of MiCPPerformance + MiCPPerformance& get_perf() {return m_perf;} + +private: + + //! \brief Return the step corrected step length if it is too long + double is_step_too_long(Eigen::VectorXd& update); + + std::shared_ptr m_program; //!< Smart pointer of a program + + MiCPSolverOptions m_options; //!< The options + MiCPPerformance m_perf; //!< The performance + + // Residuals and jacobian + Eigen::VectorXd m_residuals; //!< The residuals + Eigen::VectorXd m_phi_residuals; //!< The reformulated residuals + Eigen::VectorXd m_grad_phi; //!< The gradient of the reformulated residuals + Eigen::MatrixXd m_jacobian; //!< The jacobian + + std::vector m_max_merits; //!< Contains the m best value of the merit function + + bool m_max_taken; //!< True if the 'length' of the step was equal or bigger than the maximum step length + int m_consec_max; //!< The number of consecutive step were the step length was equal or bigger than the maximum step length + + bool m_gradient_step_taken; //!< True if the update was computed using the gradient + +}; + +} // end namespace micpsolver +} // end namespace specmicp + +// ###############// +// Implementation // +// ###############// +#include "micpsolverold.inl" + +#endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP diff --git a/src/micpsolver/micpsolverold.inl b/src/micpsolver/micpsolverold.inl new file mode 100644 index 0000000..c242f3b --- /dev/null +++ b/src/micpsolver/micpsolverold.inl @@ -0,0 +1,613 @@ +/*------------------------------------------------------- + + - Module : micpsolver + - File : micpsolver.inl + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#include "micpsolverold.hpp" // for syntaxic coloration... + +#include "estimate_cond_number.hpp" +#include "utils/log.hpp" + +#include + +//! \file micpsolver.inl implementation of the MiCP solver + +namespace specmicp { +namespace micpsolver { + +// Main algorithm +// ############## + +template +MiCPSolverReturnCode MiCPSolverOLD::solve(Eigen::VectorXd &x) +{ + int cnt = 0; + if (get_options().use_crashing) crashing(x); + else setup_residuals(x); + MiCPSolverReturnCode retcode = MiCPSolverReturnCode::NotConvergedYet; + Eigen::VectorXd update(get_neq()); + while (retcode == MiCPSolverReturnCode::NotConvergedYet) + { + DEBUG << "Iteration : " << cnt; + SPAM << "Solution : \n" << x; + m_program->hook_start_iteration(x, m_phi_residuals.norm()); + setup_residuals(x); + SPAM << "Residuals : \n ----- \n" << m_phi_residuals << "\n ----- \n"; + retcode = check_convergence(cnt, update, x); + get_perf().return_code = retcode; + if (retcode != MiCPSolverReturnCode::NotConvergedYet) break; + ++cnt; + m_max_taken = false; + setup_jacobian(x); + if(get_options().use_scaling) + search_direction_calculation(update); + else + search_direction_calculation_no_scaling(update); + reformulate_result(get_neq(), get_neq_free(), + x, m_residuals, + m_grad_phi, update); + int termcode = linesearch(update, x); + DEBUG << "Return LineSearch : " << termcode; + projection(x); + get_perf().nb_iterations = cnt; + } + return retcode; +} + +template +MiCPSolverReturnCode MiCPSolverOLD::check_convergence(int nb_iterations, + Eigen::VectorXd& update, + Eigen::VectorXd& solution) +{ + MiCPSolverReturnCode termcode = MiCPSolverReturnCode::NotConvergedYet; + const double norm_residuals = m_phi_residuals.lpNorm(); + if (norm_residuals < get_options().fvectol) + { + termcode = MiCPSolverReturnCode::ResidualMinimized; + } + else if (nb_iterations >0 and norm_update(update, solution) < get_options().steptol) + { + if (norm_residuals > get_options().threshold_stationary_point) + { + ERROR << "Stationary point detected !"; + termcode = MiCPSolverReturnCode::StationaryPoint; + } + WARNING << "Error is minimized - may indicate a stationnary point"; + termcode = MiCPSolverReturnCode::ErrorMinimized; + } + else if (nb_iterations > get_options().max_iter) + { + ERROR << "Maximum number of iteration reached (" << get_options().max_iter << ")"; + termcode = MiCPSolverReturnCode::MaxIterations; + } + else if (m_max_taken) + { + ++m_consec_max; + if (m_consec_max == get_options().maxiter_maxstep) { + ERROR << "Divergence detected - Maximum step length taken two many times"; + termcode = MiCPSolverReturnCode::MaxStepTakenTooManyTimes; + } + } + else + { + m_consec_max = 0; + } + return termcode; +} + + +template +MiCPSolverReturnCode MiCPSolverOLD::search_direction_calculation(Eigen::VectorXd& update) +{ + Eigen::VectorXd rscaler(Eigen::VectorXd::Ones(m_jacobian.cols())); + Eigen::VectorXd cscaler(Eigen::VectorXd::Ones(m_jacobian.rows())); + scaling_jacobian(m_jacobian, m_phi_residuals, rscaler, cscaler); + m_jacobian = rscaler.asDiagonal() * (m_jacobian) * cscaler.asDiagonal(); + Eigen::ColPivHouseholderQR solver; + m_gradient_step_taken = false; + int m; + for (m=0; m()); + if (cond > get_options().condition_limit) + { + m_gradient_step_taken = true; + m_jacobian += rscaler.asDiagonal() * ( + lambda*Eigen::MatrixXd::Identity(m_jacobian.rows(),m_jacobian.cols())) * cscaler.asDiagonal(); + continue; + } + update = solver.solve(-rscaler.cwiseProduct(m_phi_residuals + m*lambda*m_grad_phi)); + update = cscaler.cwiseProduct(update); + double descent_cond = m_grad_phi.dot(update); + double norm_grad = m_grad_phi.norm(); + double norm_update = update.norm(); + + if ( (descent_cond <= -get_options().factor_descent_condition*std::min(std::pow(norm_update,2),std::pow(norm_update,3))) + and (descent_cond <= -get_options().factor_descent_condition*std::min(std::pow(norm_grad,2),std::pow(norm_grad,3))) + ) + break; // we have a solution ! + m_gradient_step_taken = true; + m_jacobian += rscaler.asDiagonal() * ( lambda* + Eigen::MatrixXd::Identity(m_jacobian.rows(),m_jacobian.cols()) + ) * cscaler.asDiagonal(); + } + DEBUG << "Gradient step : m = " << m; + if (m == get_options().max_factorization_step) { + INFO << "Full gradient step taken !"; + update = -m_grad_phi; + } + return MiCPSolverReturnCode::NotConvergedYet; +} + +template +MiCPSolverReturnCode MiCPSolverOLD::search_direction_calculation_no_scaling(Eigen::VectorXd& update) +{ + DEBUG << "Solving linear system"; + Eigen::ColPivHouseholderQR solver; + m_gradient_step_taken = false; + int m; + for (m=0; m()); + if (cond > get_options().condition_limit) + { + continue; + } + update = solver.solve(-(m_phi_residuals + m*lambda*m_grad_phi)); + double descent_cond = m_grad_phi.dot(update); + double norm_grad = m_grad_phi.norm(); + double norm_update = update.norm(); + + if ( (descent_cond <= -get_options().factor_descent_condition*std::min(std::pow(norm_update,2),std::pow(norm_update,3))) + and (descent_cond <= -get_options().factor_descent_condition*std::min(std::pow(norm_grad,2),std::pow(norm_grad,3))) + ) + break; // we have a solution ! + m_gradient_step_taken = true; + m_jacobian += lambda*Eigen::MatrixXd::Identity(m_jacobian.rows(),m_jacobian.cols()); + } + DEBUG << "Gradient step : m = " << m; + if (m ==4) { + INFO << "Full gradient step taken !"; + update = -m_grad_phi; + } + return MiCPSolverReturnCode::NotConvergedYet; +} + +template +void MiCPSolverOLD::crashing(Eigen::VectorXd &x) +{ + DEBUG << "Crashing "; + const double beta = 0.5; + const double sigma = 1e-5; + int cnt = 0; + while (cnt < 10) + { + setup_residuals(x); + setup_jacobian(x); + m_grad_phi = m_jacobian.transpose()*m_phi_residuals; + Eigen::VectorXd xp(get_neq()); + int l=0; + const int maxl = 10; + while (l +void MiCPSolverOLD::reformulate_residuals(const Eigen::VectorXd& x, + const Eigen::VectorXd& r, + Eigen::VectorXd& r_phi) +{ + r_phi.resizeLike(r); + r_phi.block(0, 0, get_neq_free(), 1) = r.block(0, 0, get_neq_free(), 1); + for (int i = get_neq_free(); i +void MiCPSolverOLD::reformulate_residuals_inplace(const Eigen::VectorXd& x, + Eigen::VectorXd& r) +{ + for (int i = get_neq_free(); i +void MiCPSolverOLD:: scaling_jacobian( + const Eigen::MatrixXd& jacobian, + const Eigen::VectorXd& r_phi, + Eigen::VectorXd& rscaler, + Eigen::VectorXd& cscaler) +{ + for (int i=0; i +int MiCPSolverOLD::linesearch(Eigen::VectorXd& p, Eigen::VectorXd& x) +{ + // Reference Algo A6.3.1 : Dennis and Schnabel (1983) + DEBUG << "Linesearch"; + Eigen::VectorXd xp(get_neq()); + Eigen::VectorXd new_res(get_neq()); + double fcp; + + m_max_taken = false; + int retcode = 2; + const double alpha = 1e-6; + double newtlen = is_step_too_long(p); + double init_slope = m_grad_phi.dot(p); + double rellength = std::abs(p(0)); + for (int i=1; imax_lambda(x, p); + double lambda_prev = lambda; + + // non monotone linesearch + // ----------------------- + + double merit_value = 0.5*m_phi_residuals.squaredNorm(); + // new residual + xp = x + lambda*p; + compute_residuals(xp, new_res); + reformulate_residuals_inplace(xp, new_res); + fcp = 0.5*new_res.squaredNorm(); + + // Skip linesearch if enough progress is done + if (fcp < get_options().coeff_accept_newton_step *merit_value) + { + if (m_max_merits.size() > 0) m_max_merits[m_max_merits.size()-1] = merit_value; + else m_max_merits.push_back(merit_value); + x = xp; + return 0; + } + + //std::cout << "Merit value : " << merit_value << std::endl; + double mmax = merit_value; + if (m_max_merits.size() > 0) + { + mmax = m_max_merits[m_max_merits.size()-1]; + } + if (m_max_merits.size() < 4) + { + m_max_merits.push_back(merit_value); + if (merit_value < mmax) merit_value = (3*merit_value + mmax)/4; + } + else if (merit_value < mmax) + { + m_max_merits[3] = merit_value; + merit_value = mmax; + } + if (m_gradient_step_taken) + { + merit_value *= 100; + } + //std::cout << "Merit value used : " << merit_value << std::endl; + double fc = merit_value; + double fcp_prev; + int cnt = 0; + do + { + fcp = 0.5*new_res.squaredNorm(); + //std::cout << "fcp : " << fcp << "\n fc+alin : " << fc+alpha*lambda*init_slope << " # fc : " << fc << std::endl; + if (fcp <= fc - std::min(-alpha*lambda*init_slope,(1-alpha)*fc)) //pg760 Fachinei2003 + { + retcode = 0; + if (lambda ==1 and (newtlen > 0.99 * get_options().maxstep)) { + m_max_taken = true; + } + break; + } + else if (lambda < minlambda) + { + retcode = 1; + break; + } + else + { + double lambdatmp; + if (cnt == 0) { // only a quadratic at the first + lambdatmp = - init_slope / (2*(fcp - fc -init_slope)); + } + else + { + const double factor = 1 /(lambda - lambda_prev); + const double x1 = fcp - fc - lambda*init_slope; + const double x2 = fcp_prev - fc - lambda_prev*init_slope; + + const double a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev)); + const double b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev)); + + if (a == 0) + { // cubic interpolation is in fact a quadratic interpolation + lambdatmp = - init_slope/(2*b); + } + else + { + const double disc = b*b-3*a*init_slope; + lambdatmp = (-b+std::sqrt(disc))/(3*a); + } + if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda; + } + lambda_prev = lambda; + fcp_prev = fcp; + if (lambdatmp < 0.1*lambda) { + lambda = 0.1 * lambda; + } else { + lambda = lambdatmp; + } + } + xp = x + lambda*p; + compute_residuals(xp, new_res); + reformulate_residuals_inplace(xp, new_res); + ++cnt; + } while(retcode == 2 and cnt < 100); + DEBUG << "Lambda : " << lambda; + if (cnt == 100) + { + ERROR << "Too much linesearch iterations ! We stop"; + } + x = xp; + p = lambda*p; + return retcode; + +} + +// Projection of the variables onto the feasible set +template +void MiCPSolverOLD::projection(Eigen::VectorXd &x) +{ + for (int i=0; inb_complementarity_variables(); ++i) + { + if (x(i+m_program->nb_free_variables()) < get_options().projection_min_variable) + { + x(i+m_program->nb_free_variables()) = 0; + } + } +} + +template +double MiCPSolverOLD::is_step_too_long(Eigen::VectorXd& update) +{ + double steplength = update.norm(); + if (steplength > get_options().maxstep) + { + m_max_taken = true; + update = get_options().maxstep / steplength * update; + steplength = get_options().maxstep; + } + return steplength; +} + + +// ================================================= // +// // +// NCP functions and reformulation // +// // +// ================================================= // + +template <> +inline double ncp_function(double a, double b, double t) +{ + return penalized_fisher_burmeister(a, b, t); +} + +template <> +inline double ncp_function(double a, double b, double _) +{ + return std::min(a, b); +} + +template <> +inline void reformulate_jacobian_helper( + int neq, + int neq_free, + const Eigen::VectorXd& x, + const Eigen::VectorXd& r, + Eigen::MatrixXd& jacobian, + Eigen::VectorXd& _, + double t + ) +{ + // set the z vector : contains 1 for degenerate points + Eigen::VectorXd z(Eigen::VectorXd::Zero(neq)); + for (int i=neq_free; i 0) and (x(i) >0)) + { + c -= (1-lambda)*r(i); + d -= (1-lambda)*x(i); + } + + grad_fi = d*grad_fi; + grad_fi(i) += c; + } + jacobian.block(i, 0, 1, neq) = grad_fi.transpose(); + } +} + +template <> +inline void reformulate_jacobian_helper( + int neq, + int neq_free, + const Eigen::VectorXd& x, + const Eigen::VectorXd& r, + Eigen::MatrixXd& jacobian, + Eigen::VectorXd& r_phi, + double _ + ) +{ + std::vector to_keep; + to_keep.reserve(10); + Eigen::VectorXd to_remove(neq-neq_free); + for (int i=neq_free; i= r(i)) + { + to_remove(i-neq_free) = 0; + to_keep.push_back(i); + } + else + to_remove(i-neq_free) = x(i); + } + r_phi.block(0, 0, neq_free, 1) -= jacobian.block(0, neq_free, neq_free, neq-neq_free)*to_remove; + int new_i = neq_free; + for (auto it=to_keep.begin(); it!=to_keep.end(); ++it) + { + //r_phi.block(0, 0, neq_free, 1) += x(*it)*jacobian.block(0, *it, neq_free, 1); + jacobian.block(new_i, 0, 1, neq_free) = jacobian.block(*it, 0, 1, neq_free); // the bottom right corner is 0 + jacobian.block(0, new_i, neq_free, 1) = jacobian.block(0, *it, neq_free, 1); + r_phi(new_i) = r_phi(*it); + ++new_i; + + } + r_phi.conservativeResize(new_i); + jacobian.conservativeResize(new_i, new_i); + DEBUG << jacobian; +} + + +template <> +inline void reformulate_result( + int neq, + int neq_free, + Eigen::VectorXd& x, + const Eigen::VectorXd& orig_r, + Eigen::VectorXd& grad_phi, + Eigen::VectorXd& update + ) +{} + +template <> +inline void reformulate_result( + int neq, + int neq_free, + Eigen::VectorXd& x, + const Eigen::VectorXd& orig_r, + Eigen::VectorXd& grad_phi, + Eigen::VectorXd& update + ) +{ + //std::cout << " Update \n ------- \n " << update << std::endl; + int tot_to_keep = 0; + for (int i=neq_free; i= orig_r(i)) + ++tot_to_keep; + } + //std::cout << " update \n ------ \n" << update.block(neq_free, 0, tot_to_keep, 1) << std::endl; + update.conservativeResize(neq); + grad_phi.conservativeResize(neq); + int kept_i = 1; + for (int i=neq-1; i>=neq_free; --i) + { // we go backwards to avoid extra copies + //std::cout << i << " # " << x(i) << " - " << orig_r(i) << std::endl; + if (x(i) >= orig_r(i)) + { + //std::cout << i << std::endl; + update(i) = update(neq_free+(tot_to_keep-kept_i)); + //std::cout << update(i) << std::endl; + grad_phi(i) = grad_phi(neq_free+(tot_to_keep-kept_i)); + ++kept_i; + } + else + { + //x(i) = 0.0; + //update(i) = 0.0; + update(i) = -x(i); + grad_phi(i) = x(i); + } + } +} + +} // end namespace micpsolver +} // end namespace specmicp diff --git a/src/specmicp/reduced_system.cpp b/src/specmicp/reduced_system.cpp index 635cfea..edec286 100644 --- a/src/specmicp/reduced_system.cpp +++ b/src/specmicp/reduced_system.cpp @@ -1,480 +1,485 @@ /*------------------------------------------------------- - Module : specmicp - File : specmicp.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 "reduced_system.hpp" #include "utils/log.hpp" #include "equilibrium_data.hpp" #include "physics/constants.hpp" #include "physics/laws.hpp" namespace specmicp { inline double pow10(double x) { return std::pow(10.0, x); } ReducedSystem::ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc): m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(ptrthermo->nb_aqueous()), m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous()) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(true); m_secondary_conc.setZero(); m_loggamma.setZero(); } ReducedSystem::ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc, const SimulationBox &simul_box): m_simulbox(simul_box), m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(ptrthermo->nb_aqueous()), m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous()) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(true); m_secondary_conc.setZero(); m_loggamma.setZero(); } ReducedSystem::ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc, const SimulationBox& simul_box, const EquilibriumState& previous_solution): m_simulbox(simul_box), m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(previous_solution.molality_secondary()), m_loggamma(previous_solution.activity_coefficient()) { if (not m_simulbox.is_fixed_temperature() or m_simulbox.get_temperature() != units::celsius(25)) { throw std::runtime_error("Only a fixed temperature of 25°C is supported"); } number_eq(true); } void ReducedSystem::number_eq(bool water_equation) { int neq = 0; m_ideq.reserve(m_thermo->nb_component()+m_thermo->nb_mineral()); if (water_equation) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } for (int i=1;inb_component();++i) { // Remove components that does not exist if (m_tot_conc(i) == 0 and i!= 1) { INFO << "Component " << m_thermo->label_aqueous(i) << "has total concentration equal to zero, we desactivate it" << std::endl; m_nonactive_component.push_back(i); m_ideq.push_back(no_equation); continue; } else { m_ideq.push_back(neq); ++neq; } } m_nb_free_vars = neq; for (int m=0; mnb_mineral();++m) { bool can_precipitate = true; // Remove minerals that cannot precipitate for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_thermo->nu_min(m, *it) != 0) { can_precipitate = false; break; // this is not a mineral that can precipitate } } if (can_precipitate) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } } m_nb_tot_vars = neq; m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars; // aqueous species m_active_aqueous.reserve(m_thermo->nb_aqueous()); for (int j=0; jnb_aqueous(); ++j) { bool can_exist = true; for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_thermo->nu_aq(j,*it) != 0) { can_exist = false; } } m_active_aqueous.push_back(can_exist); } } double ReducedSystem::residual_water(const Eigen::VectorXd& x) const { double res = m_tot_conc(0) - mass_water(x)/molar_mass_water; for (int j=0; jnb_aqueous(); ++j) { + if (m_thermo->nu_aq(j, 0) == 0) continue; res -= mass_water(x)*m_thermo->nu_aq(j, 0)*m_secondary_conc(j); } for (int m=0; mnb_mineral(); ++m) { - if (ideq_min(m) == no_equation) continue; + if (ideq_min(m) == no_equation or m_thermo->nu_min(m, 0) == 0) continue; res -= m_thermo->nu_min(m, 0)*x(ideq_min(m)); } return res; } double ReducedSystem::residual_component(const Eigen::VectorXd& x, int i) const { const int idp = ideq_paq(i); double res = m_tot_conc(i) - mass_water(x)*pow10(x(idp)); for (int j=0; jnb_aqueous(); ++j) { if (not m_active_aqueous[j]) continue; res -= mass_water(x)*m_thermo->nu_aq(j, i)*m_secondary_conc(j); } for (int m=0; mnb_mineral(); ++m) { if (ideq_min(m) == no_equation) continue; res -= m_thermo->nu_min(m, i)*x(ideq_min(m)); } return res; } double ReducedSystem::residual_mineral(const Eigen::VectorXd& x, int m) const { double res = m_thermo->logkr_min(m); for (int i=1; inb_component(); ++i) { if (m_thermo->nu_min(m, i) != 0) res -= m_thermo->nu_min(m, i)*(x(ideq_paq(i)) + m_loggamma(i)); } return res; } void ReducedSystem::get_residuals(const Eigen::VectorXd& x, Eigen::VectorXd& residual) { set_secondary_concentration(x); if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x); for (int i=1; inb_component(); ++i) { if (ideq_paq(i) != no_equation) residual(ideq_paq(i)) = residual_component(x, i); } for (int m=0; mnb_mineral(); ++m) { if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m); } } //non-optimized Finite difference, for test only //void ReducedSystem::get_jacobian(Eigen::VectorXd& x, // Eigen::MatrixXd& jacobian) //{ // const int neq = total_variables(); // Eigen::VectorXd res(total_variables()); // Eigen::VectorXd perturbed_res(total_variables()); // get_residuals(x, res); // for (int j=0; jnb_aqueous(); ++j) { if ( m_thermo->nu_aq(j, 0) == 0 ) continue; tmp -= m_thermo->nu_aq(j, 0)*m_secondary_conc(j); } jacobian(0, 0) = tmp; for (int k=1; knb_component(); ++k) { if (ideq_paq(k) == no_equation) continue; double tmp = 0; for (int j=0; jnb_aqueous(); ++j) { tmp -= m_thermo->nu_aq(j,0)*m_thermo->nu_aq(j,k)*m_secondary_conc(j); } jacobian(0, ideq_paq(k)) = x(ideq_w())*log10 * tmp; } for (int m=0; mnb_mineral(); ++m) { if (ideq_min(m) == no_equation) continue; jacobian(0, ideq_min(m)) = -m_thermo->nu_min(m, 0); } } // aqueous component for (int i=1; inb_component(); ++i) { if (ideq_paq(i) == no_equation) continue; const int idp = ideq_paq(i); // aqueous components for (int k=1; knb_component(); ++k) { if (ideq_paq(k) == no_equation) continue; double tmp_iip = 0; if (k == i) tmp_iip -= pow10(x(ideq_paq(i)))*log10; for (int j=0; jnb_aqueous(); ++j) { tmp_iip -= log10*m_thermo->nu_aq(j, i)*m_thermo->nu_aq(j, k)*m_secondary_conc(j); } jacobian(idp, ideq_paq(k)) = mass_water(x)*tmp_iip; } // minerals for (int m=0; mnb_mineral(); ++m) { if (ideq_min(m) == no_equation) continue; jacobian(idp, ideq_min(m)) = - m_thermo->nu_min(m, i); } if (ideq_w() != no_equation) // Water { double tmp_iw = -std::pow(10.0, x(idp)); for (int j=0; jnb_aqueous(); ++j) { if ( m_thermo->nu_aq(j, i) == 0 ) continue; tmp_iw -= m_thermo->nu_aq(j, i)*m_secondary_conc(j); } jacobian(idp, ideq_w()) = tmp_iw; } } // mineral equilibrium for (int m=0; mnb_mineral(); ++m) { const int idm = ideq_min(m); if (idm == no_equation) continue; for (int i=1; inb_component(); ++i) { if (ideq_paq(i) == no_equation) continue; jacobian(idm, ideq_paq(i)) = -m_thermo->nu_min(m, i); } } } void ReducedSystem::set_secondary_concentration(const Eigen::VectorXd& x) { for (int j=0; jnb_aqueous(); ++j) { if (m_active_aqueous[j] == false) { m_secondary_conc(j) = 0; continue; } double conc = -m_thermo->logkr_aq(j) - m_loggamma(j+m_thermo->nb_component()); for (int k=1; knb_component(); ++k) { if (m_thermo->nu_aq(j, k) == 0) continue; conc += m_thermo->nu_aq(j, k)*(x(ideq_paq(k)) + m_loggamma(k)); } m_secondary_conc(j) = pow10(conc); } } void ReducedSystem::set_ionic_strength(const Eigen::VectorXd& x) { double ionic = 0; for (int i=1; inb_component(); ++i) { if (ideq_paq(i) == no_equation or m_thermo->charge_paq(i) == 0) continue; ionic += pow10(x(ideq_paq(i)))*std::pow(m_thermo->charge_paq(i),2); } for (int j=0; jnb_aqueous(); ++j) { if (not m_active_aqueous[j] or m_thermo->charge_saq(j) == 0) continue; ionic += m_secondary_conc(j)*std::pow(m_thermo->charge_saq(j),2); } m_ionic_strength = ionic; } void ReducedSystem::compute_log_gamma(const Eigen::VectorXd& x) { set_ionic_strength(x); const double sqrti = std::sqrt(m_ionic_strength); for (int i=0; inb_component(); ++i) { if (ideq_paq(i) == no_equation) { m_loggamma(i) = 0; continue; } m_loggamma(i) = laws::extended_debye_huckel(m_ionic_strength, sqrti, m_thermo->charge(i), m_thermo->ao_debye(i), m_thermo->bdot_debye(i)); } for (int j=0; jnb_aqueous(); ++j) { if (not m_active_aqueous[j]) { m_loggamma(m_thermo->nb_component() + j) = 0; continue; } m_loggamma(m_thermo->nb_component() + j) = laws::extended_debye_huckel( m_ionic_strength, sqrti, m_thermo->charge_saq(j), m_thermo->ao_debye_saq(j), m_thermo->bdot_debye_saq(j)); } } bool ReducedSystem::hook_start_iteration(const Eigen::VectorXd& x, double norm_residual) { - //compute_log_gamma(x); not_in_linesearch = true; double previous_norm = m_loggamma.norm(); bool may_have_converged = false; if (norm_residual < nb_free_variables()) { // Use fix point iterations for non-ideality for (int i=0; i<6; ++i) { set_secondary_concentration(x); compute_log_gamma(x); - //std::cout << " i : " << i << std::endl; if (std::abs(previous_norm - m_loggamma.norm()) < 1e-6) {may_have_converged = true; break;} previous_norm = m_loggamma.norm(); } } return may_have_converged; } double ReducedSystem::max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update) { if (ideq_w() != no_equation) { - return 1.0/std::max(1.0, -update(0)/(0.25*x(0))); + return 1.0/std::max(1.0, -update(0)/(0.9*x(0))); } else { return 1.0; } } -void ReducedSystem::reasonable_starting_guess(Eigen::VectorXd& x) +void ReducedSystem::reasonable_starting_guess(Eigen::VectorXd& x, bool for_min) { x.resize(m_thermo->nb_component()+ m_thermo->nb_mineral()); - x(0) = m_tot_conc(0)*specmicp::molar_mass_water; + x(0) = 2*m_tot_conc(0)*specmicp::molar_mass_water; for (int i=1; inb_component(); ++i) { - x(i) = -2;//std::log10(10*m_tot_conc(i)); + x(i) = -4.0; } - x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero(); + if (for_min) + x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero(); + else + x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero(); + + m_loggamma.setZero(); + m_secondary_conc.setZero(); } void ReducedSystem::reasonable_restarting_guess(Eigen::VectorXd& x) { x(0) = m_tot_conc(0)*specmicp::molar_mass_water; for (int i=m_thermo->nb_component(); inb_component()+m_thermo->nb_mineral(); ++i) { - x(i) = 0.0; + x(i) = -4.0; } m_loggamma.setZero(); m_secondary_conc.setZero(); } EquilibriumState ReducedSystem::get_solution(const Eigen::VectorXd &xtot, const Eigen::VectorXd& x) { double previous_norm = m_loggamma.norm(); set_secondary_concentration(x); compute_log_gamma(x); if (std::abs(previous_norm - m_loggamma.norm()) > 1e-6) { WARNING << "Activity coefficient have not converged !" << std::endl << "output can not be trusted\n Difference "+std::to_string(std::abs(previous_norm - m_loggamma.norm())); } return EquilibriumState(xtot, m_secondary_conc, m_loggamma, m_ionic_strength, m_thermo->get_raw_database()); } } // end namespace specmicp diff --git a/src/specmicp/reduced_system.hpp b/src/specmicp/reduced_system.hpp index 5b62d9e..5b63d1f 100644 --- a/src/specmicp/reduced_system.hpp +++ b/src/specmicp/reduced_system.hpp @@ -1,203 +1,203 @@ /*------------------------------------------------------- - Module : specmicp - File : specmicp.hpp - 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. ---------------------------------------------------------*/ #ifndef SPECMICP_REDUCEDSYSTEM_HPP #define SPECMICP_REDUCEDSYSTEM_HPP //! \file reduced_system.hpp The MiCP program to solve speciation #include #include "thermodata.hpp" #include "kineticsdata.hpp" #include "micpsolver/micpprog.hpp" #include "simulation_box.hpp" //! \namespace specmicp main namespace for the SpecMiCP solver namespace specmicp { class EquilibriumState; //! \brief The equilibrium speciation problem //! //! Represent the equilibrium speciation problem as a Mixed Complementarity program //! It is called reduced becaused it does not solve explicitely for the secondary species. //! //! Equations //! ========== //! //! Main equations //! -------------- //! //! - Water conservation equation : \f[ T_w = m_w \left( \frac{1}{M_w} + \sum_{j=1}^{N_r} \nu_{wj} m_j \right) + \sum_{l=1}^{N_{min}} \nu_{lw} n_l \f] //! - Component conservation equations : \f[ T_i = m_w \left( m_i + \sum_{j=1}^{N_r} \nu_{ij} m_j \right) + \sum_{l=1}^{N_{min}} \nu_{li} n_l \f] //! - Mineral equilibrium conditions : \f[ n_l \geq 0, \; 1 - \mathrm{SI}_l \geq 0 , \; \mathrm{and} \; n_l ( 1 - \mathrm{SI}_l ) = 0 \f] //! //! Secondary equations //! ------------------- //! //! - Secondary aqueous species equilibrium conditions : \f[ m_j = \frac{\prod_{i=2}^{N_c} (\gamma_i m_i)^{\nu_{ji}}}{K_j \gamma_j} \f] //! - Ionic Strength : \f[ I = \sum_{k=2}^{N_c+N_r} z_k^2 m_k \f] //! - Activity coefficients : \f[ \log \gamma_k = \frac{-A z_k^2 \sqrt{I}}{1 + B \dot{a}_k \sqrt{I}} + \dot{b}_k I \f] //! class ReducedSystem : public micpsolver::MiCPProg { public: const int no_equation = -1; //! The corresponding variables is not a degree of freedom (i.e. is fixed) //! \brief Create a reduced system //! //! \param ptrthermo Shared pointer to the thermodynamic data //! \param totconc Vector containing the total concentrations (in mol) for each components ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc); //! \brief Create a reduced system //! //! \param ptrthermo Shared pointer to the thermodynamic data //! \param totconc Vector containing the total concentrations (in mol) for each components ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc, const SimulationBox& simul_box); ReducedSystem(std::shared_ptr ptrthermo, const Eigen::VectorXd& totconc, const SimulationBox& simul_box, const EquilibriumState& previous_solution); // Variables // ========== //! \brief Return the total number of variables int total_variables() const {return m_nb_tot_vars;} //! \brief Return the number of 'free' variables (not subject to complementarity conditions) int nb_free_variables() const {return m_nb_free_vars;} //! \brief Return the number of variables subjected to the complementarity conditions int nb_complementarity_variables() const {return m_nb_compl_vars;} // The linear system // ================== //! \brief Return the residual for the water conservation equation double residual_water(const Eigen::VectorXd& x) const; //! \brief Return the residual for the conservation equation of component (i) //! //! For (i==0), water use the method : 'residual_water' double residual_component(const Eigen::VectorXd& x, int i) const; //! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$ double residual_mineral(const Eigen::VectorXd& x, int m) const; //! \brief Return the residuals void get_residuals(const Eigen::VectorXd& x, Eigen::VectorXd& residual); //! \brief Return the jacobian void get_jacobian(Eigen::VectorXd &x, Eigen::MatrixXd& jacobian); //! \brief Return the equation id int ideq(int i) const {return m_ideq[i];} //! \brief Return the equation number for conservation of water int ideq_w() const {return m_ideq[0];} //! \brief Return the equation number of component 'i' int ideq_paq(int i) const { assert(i < m_thermo->nb_component()); return m_ideq[i]; } //! \brief Return the equation number of mineral 'm' int ideq_min(int m) const { assert(m < m_thermo->nb_mineral()); return m_ideq[m+m_thermo->nb_component()]; } //! \brief Return the mass of water double mass_water(const Eigen::VectorXd& x) const { if (ideq_w() != no_equation) return x(ideq_w()); else return 1.0; } //! \brief Compute the ionic strength void set_ionic_strength(const Eigen::VectorXd& x); //! \brief Compute the activity coefficients void compute_log_gamma(const Eigen::VectorXd& x); //! \brief Compute the activity model, called by the solver at the beginning of an iteration bool hook_start_iteration(const Eigen::VectorXd& x, double norm_residual); //! \brief Return a scale factor to avoid negative mass during Newton's iteration double max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update); //! \brief Return the component that does not exist in the solution (inactive dof) const std::vector& get_non_active_component() {return m_nonactive_component;} //! \brief A reasonable (well... maybe...) starting guess - void reasonable_starting_guess(Eigen::VectorXd& x); + void reasonable_starting_guess(Eigen::VectorXd& x, bool for_min=false); //! \brief A reasonable (maybe...) restarting guess void reasonable_restarting_guess(Eigen::VectorXd& x); //! Return the equilibrium state of the system, the 'Solution' of the speciation problem EquilibriumState get_solution(const Eigen::VectorXd& xtot, const Eigen::VectorXd& x); private: // For the normal algorithm void set_secondary_concentration(const Eigen::VectorXd& x); // For the tot aq_conc pb //void unsafe_set_secondary_concentration(const Eigen::VectorXd& x); void number_eq(bool water_equation); SimulationBox m_simulbox; std::shared_ptr m_thermo; Eigen::VectorXd m_tot_conc; int m_nb_tot_vars; int m_nb_free_vars; int m_nb_compl_vars; std::vector m_ideq; Eigen::VectorXd m_secondary_conc; Eigen::VectorXd m_gas_pressure; double m_ionic_strength; Eigen::VectorXd m_loggamma; bool not_in_linesearch; std::vector m_nonactive_component; std::vector m_active_aqueous; std::vector m_active_gas; }; } // end namespace specmicp #endif // SPECMICP_REDUCEDSYSTEM_HPP diff --git a/src/specmicp/reduced_system_solver.cpp b/src/specmicp/reduced_system_solver.cpp index 7a575fc..21408dd 100644 --- a/src/specmicp/reduced_system_solver.cpp +++ b/src/specmicp/reduced_system_solver.cpp @@ -1,218 +1,226 @@ /*------------------------------------------------------- - 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 "reduced_system_solver.hpp" #include "micpsolver/micpsolver.hpp" -#include "micpsolver/micpsolver_min.hpp" + +#include "micpsolver/micpsolverold.hpp" +//#include "micpsolver/micpsolver_min.hpp" #include "equilibrium_data.hpp" #include "simulation_box.hpp" #include namespace specmicp { ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr data, const Eigen::VectorXd& totconc, const SimulationBox& simulbox): m_options(), m_data(data), m_system(std::make_shared(std::make_shared(data), totconc, simulbox)), m_var(Eigen::VectorXd::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr data, const Eigen::VectorXd& totconc, const SimulationBox& simulbox, const EquilibriumState& previous_solution): m_options(), m_data(data), m_system(std::make_shared(std::make_shared(data), totconc, simulbox, previous_solution)), m_var(Eigen::VectorXd::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr data, const Eigen::VectorXd& totconc): m_options(), m_data(data), m_system(std::make_shared(std::make_shared(data), totconc)), m_var(Eigen::VectorXd::Zero(data->nb_component+data->nb_mineral)) {} ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr data, const Eigen::VectorXd& totconc, const EquilibriumState& previous_solution): m_options(), m_data(data), m_system(std::make_shared(std::make_shared(data), totconc, SimulationBox(), previous_solution)), m_var(Eigen::VectorXd::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(Eigen::VectorXd &x, bool init) { - if (init) m_system->reasonable_starting_guess(x); + 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; 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 Eigen::VectorXd &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 m_var = x; // direct copy for (int i=0; inb_component; ++i) { 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; } assert(new_i == m_data->nb_component - non_active_component.size()); for (int m=0; mnb_mineral; ++m) { 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); } } 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::MiCPSolverMin solver(m_system); + 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(Eigen::VectorXd &x) { const std::vector& non_active_component = m_system->get_non_active_component(); if (non_active_component.size() == 0) { x = m_var; //direct copy // at that point we should have the correct solution } else { unsigned int new_i = 0; for (int i=0; inb_component; ++i) { 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; } assert(new_i == m_data->nb_component - non_active_component.size()); for (int m=0; mnb_mineral; ++m) { 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; } } } } } // end namespace specmicp