diff --git a/src/specmicp/reduced_system.cpp b/src/specmicp/reduced_system.cpp index 7825b23..ebc8b69 100644 --- a/src/specmicp/reduced_system.cpp +++ b/src/specmicp/reduced_system.cpp @@ -1,540 +1,541 @@ /*------------------------------------------------------- - Module : specmicp - File : specmicp.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #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, Eigen::VectorXd& totconc): m_include_gas_phase(false), m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(ptrthermo->nb_aqueous()), m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous()) { number_eq(true); m_secondary_conc.setZero(); m_loggamma.setZero(); 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"); } } ReducedSystem::ReducedSystem(std::shared_ptr ptrthermo, Eigen::VectorXd& totconc, SimulationBox simul_box): m_simulbox(simul_box), m_include_gas_phase(false), m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(ptrthermo->nb_aqueous()), m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous()) { number_eq(true); m_secondary_conc.setZero(); m_loggamma.setZero(); if (m_simulbox.is_fixed_volume() and m_thermo->nb_gas() != 0) { m_include_gas_phase = true; } 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"); } } ReducedSystem::ReducedSystem(std::shared_ptr ptrthermo, Eigen::VectorXd& totconc, SimulationBox simul_box, const EquilibriumState& previous_solution): m_simulbox(simul_box), m_include_gas_phase(false), m_thermo(ptrthermo), m_tot_conc(totconc), m_secondary_conc(previous_solution.molality_secondary()), m_loggamma(previous_solution.activity_coefficient()) { number_eq(true); if (m_simulbox.is_fixed_volume() and m_thermo->nb_gas() != 0) { m_include_gas_phase = true; } 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"); } } 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); } // gas if (is_gas_phase_computed()) { m_active_gas.reserve(m_thermo->nb_gas()); for (int g=0; gnb_gas(); ++g) { bool can_exist = true; for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_thermo->nu_gas(g,*it) != 0) { can_exist = false; } } m_active_gas.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) { 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; res -= m_thermo->nu_min(m, 0)*x(ideq_min(m)); } if (is_gas_phase_computed()) { for (int g=0; gnb_gas(); ++g) { res -= m_thermo->nu_gas(g, 0)*laws::mole_perfect_gas( m_gas_pressure(g),m_tmp_gas_volume , m_simulbox.get_temperature()); } } 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)); } if (is_gas_phase_computed()) { for (int g=0; gnb_gas(); ++g) { res -= m_thermo->nu_gas(g, i)*laws::mole_perfect_gas( m_gas_pressure(g), m_tmp_gas_volume, m_simulbox.get_temperature()); } } 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); } } //std::cout << jacobian << std::endl; } 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); } if (is_gas_phase_computed()) { for (int g=0; gnb_gas(); ++g) { if (m_active_gas[g] == false) { m_gas_pressure(g) = 0; continue; } double logpres = -m_thermo->logkr_gas(g); for (int k=1; knb_component(); ++k) { if (m_thermo->nu_gas(g, k) == 0) continue; logpres += m_thermo->nu_gas(g, k)*(x(ideq_paq(k)) +m_loggamma(k)); } m_gas_pressure(g) = pow10(logpres); } } } + 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(); } } //std::cout << may_have_converged << std::endl; return may_have_converged; } double ReducedSystem::max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update) { //return 1.0/std::max(1.0,(-update.block(0,0,5,1).cwiseQuotient(0.95*x.block(0,0,5,1))).maxCoeff()); if (ideq_w() != no_equation) { // std::cout << "max lambda : " << 1.0/std::max(1.0, -update(0)/(0.95*x(0))) << std::endl; return 1.0/std::max(1.0, -update(0)/(0.25*x(0))); } else { return 1.0; } } void ReducedSystem::reasonable_starting_guess(Eigen::VectorXd& x) { x.resize(total_variables()); x(0) = m_tot_conc(0)*specmicp::molar_mass_water; for (int i=1; inb_component(); ++i) { x(i) = std::log10(0.1*m_tot_conc(i)/x(0)); } x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero(); } void ReducedSystem::reasonable_restarting_guess(Eigen::VectorXd& x) { //x.conservativeResize(total_variables()); x(0) = m_tot_conc(0)*specmicp::molar_mass_water; //for (int i=1; inb_component(); ++i) //{ // x(i) = std::log10(0.1*m_tot_conc(i)/x(0)); //} //x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero(); for (int i=m_thermo->nb_component(); inb_component()+m_thermo->nb_mineral(); ++i) { x(i) = 0.0; } // reset secondary variable m_loggamma.setZero(); m_secondary_conc.setZero(); } -EquilibriumState ReducedSystem::get_solution(const Eigen::VectorXd& x) +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 - output can not be trusted\n Difference "+std::to_string(std::abs(previous_norm - m_loggamma.norm())); } - return EquilibriumState(x, m_secondary_conc, m_loggamma, m_ionic_strength, m_thermo->get_raw_database()); + return EquilibriumState(xtot, m_secondary_conc, m_loggamma, m_ionic_strength, m_thermo->get_raw_database()); } double ReducedSystem::condensed_phase_volume(const Eigen::VectorXd& x) const { double volume = x(ideq_w())/constants::water_density_25; for (int m=0; mnb_mineral(); ++m) { volume += x(ideq_min(m))*m_thermo->get_raw_database()->molar_volume_mineral(m); } return volume; } double ReducedSystem::gasphase_volume(const Eigen::VectorXd& x) const { assert(m_simulbox.is_fixed_volume()); return (m_simulbox.get_volume() - condensed_phase_volume(x)); } } // end namespace specmicp diff --git a/src/specmicp/reduced_system.hpp b/src/specmicp/reduced_system.hpp index ae619f8..75c6331 100644 --- a/src/specmicp/reduced_system.hpp +++ b/src/specmicp/reduced_system.hpp @@ -1,188 +1,188 @@ /*------------------------------------------------------- - Module : specmicp - File : specmicp.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #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, 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, Eigen::VectorXd& totconc, SimulationBox simul_box); ReducedSystem(std::shared_ptr ptrthermo, Eigen::VectorXd& totconc, 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); //! \brief A reasonable (maybe...) restarting guess void reasonable_restarting_guess(Eigen::VectorXd& x); - EquilibriumState get_solution(const Eigen::VectorXd& x); + EquilibriumState get_solution(const Eigen::VectorXd& xtot, const Eigen::VectorXd& x); //! Return true if the computation include the gas phase bool is_gas_phase_computed() const {return m_include_gas_phase;} 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); double condensed_phase_volume(const Eigen::VectorXd& x) const; double gasphase_volume(const Eigen::VectorXd& x) const; SimulationBox m_simulbox; bool m_include_gas_phase; 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; double m_tmp_gas_volume; }; } // end namespace specmicp #endif // SPECMICP_REDUCEDSYSTEM_HPP diff --git a/src/specmicp/reduced_system_solver.cpp b/src/specmicp/reduced_system_solver.cpp index 923b3e0..3dbdebd 100644 --- a/src/specmicp/reduced_system_solver.cpp +++ b/src/specmicp/reduced_system_solver.cpp @@ -1,174 +1,175 @@ /*------------------------------------------------------- - Module : specmicp - File : reduced_system_solver - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "reduced_system_solver.hpp" #include "micpsolver/micpsolver.hpp" #include "micpsolver/micpsolver_min.hpp" #include "equilibrium_data.hpp" #include "simulation_box.hpp" #include namespace specmicp { ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr data, 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, 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) { - return m_system->get_solution(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); 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) { m_var = x; 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) { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_performance(); } else { //micpsolver::MiCPSolver solver(m_system); micpsolver::MiCPSolverMin 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; } 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; // FIXME 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