Homec4science
Diffusion SpecMiCP / ReactMiCP (fixpoint) Browse

specmicp/src/specmicp/reaction_path.cpp
c69b6bbe42b8fixpoint

Last Change
Raw File

reaction_path.cpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : reaction_path.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "reaction_path.hpp"
#include "thermodata.hpp"
#include "micpsolver/micpsolver.hpp"
#include "utils/log.hpp"
#include <iostream>
namespace specmicp {
void ReactionPathDriver::read_database()
{
m_database.parse_database(m_model->database_path);
m_database.make_canonical();
m_data = m_database.get_database();
}
int ReactionPathDriver::label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
{
auto it = std::find(list_labels.begin(), list_labels.end(), label);
if (it == list_labels.end())
{
throw std::invalid_argument("Unknown species : " + label);
}
return it - list_labels.begin();
}
void ReactionPathDriver::dissolve_to_components()
{
m_tot_conc = Eigen::VectorXd::Zero(m_data->nb_component);
m_tot_conc_increase = Eigen::VectorXd::Zero(m_data->nb_component);
// everything is dissolved ...
for (auto it=m_model->amount_components.begin(); it!=m_model->amount_components.end(); ++it)
{
const int id = label_to_id(it->first, m_data->labels_basis);
m_tot_conc(id) = it->second.first;
m_tot_conc_increase(id) = it->second.second;
}
// Fixme : add contribution for minerals and aqueous species
std::vector<int> to_remove;
to_remove.reserve(10);
int new_i = 0;
for (int i=0; i<m_data->nb_component; ++i)
{
if (m_tot_conc(i) == 0 and m_tot_conc_increase(i) == 0)
{
to_remove.push_back(i);
continue;
}
m_tot_conc(new_i) = m_tot_conc(i);
m_tot_conc_increase(new_i) = m_tot_conc_increase(i);
++new_i;
}
m_tot_conc.conservativeResize(new_i);
m_tot_conc_increase.conservativeResize(new_i);
if (to_remove.size() > 0) m_database.remove_components(to_remove);
}
micpsolver::MiCPPerformance ReactionPathDriver::one_step(Eigen::VectorXd& x, bool init)
{
m_current_solver = ReducedSystemSolver(m_data, m_tot_conc);
m_current_solver.get_options() = m_solver_options;
micpsolver::MiCPPerformance perf = m_current_solver.solve(x, init);
if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
{
ERROR << "Failed to solve the system !";
}
m_tot_conc += m_tot_conc_increase;
return perf;
}
void ReactionPathDriver::run_all()
{
Eigen::VectorXd x;
// first step
one_step(x, true);
for (int i=1; i<m_model->nb_step; ++i)
{
one_step(x);
save_data(x);
}
}
void ReactionPathDriver::save_data(Eigen::VectorXd& x)
{
}
} // end namespace specmicp
c4science · Help