Homec4science
Diffusion SpecMiCP / ReactMiCP (specmicp_bin) Browse

specmicp/src/specmicp/reduced_system_solver.cpp
e4061f7461edspecmicp_bin

Last Change
Raw File

reduced_system_solver.cpp
View Options

/*-------------------------------------------------------
- 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 <iostream>
namespace specmicp {
ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr<database::DataContainer> data,
Eigen::VectorXd& totconc):
m_options(),
m_data(data),
m_system(std::make_shared<ReducedSystem>(std::make_shared<ThermoData>(data), totconc)),
m_var(Eigen::VectorXd::Zero(data->nb_component+data->nb_mineral))
{}
ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr<database::DataContainer> data,
Eigen::VectorXd& totconc,
const EquilibriumState& previous_solution):
m_options(),
m_data(data),
m_system(std::make_shared<ReducedSystem>(std::make_shared<ThermoData>(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);
}
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<int>& non_active_component = m_system->get_non_active_component();
if (non_active_component.size() == 0)
{
m_var = x;
for (int i=0; i<m_var.rows(); ++i)
{
if (m_var(i) == -HUGE_VAL)
{
m_var(i) = -4;
}
}
}
else
{
unsigned int new_i = 0;
for (int i=0; i<m_data->nb_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; m<m_data->nb_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<ReducedSystem> solver(m_system);
solver.set_options(get_options().solver_options);
solver.solve(m_var);
return solver.get_performance();
}
else
{
//micpsolver::MiCPSolver<ReducedSystem,micpsolver::NCPfunction::min> solver(m_system);
micpsolver::MiCPSolverMin<ReducedSystem> 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<int>& 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; i<m_data->nb_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; m<m_data->nb_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
c4science · Help