Homec4science
Diffusion SpecMiCP / ReactMiCP (gasphase) Browse

specmicp/src/specmicp/reaction_path.cpp
881a9cd5b582gasphase

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();
}
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_aqueous.begin(); it!=m_model->amount_aqueous.end(); ++it)
{
try
{
int id = database::Database::safe_label_to_id(it->first, m_data->labels_basis);
m_tot_conc(id) += it->second.first;
m_tot_conc_increase(id) += it->second.second;
}
catch(std::invalid_argument& e)
{
int id = database::Database::safe_label_to_id(it->first, m_data->labels_aqueous);
for (int idc=0; idc<m_data->nb_component; ++idc)
{
m_tot_conc(idc) += m_data->nu_aqueous(id, idc)*it->second.first;
m_tot_conc_increase(id) += m_data->nu_aqueous(id, idc)*it->second.second;
}
}
}
for (auto it=m_model->amount_minerals.begin(); it!=m_model->amount_minerals.end(); ++it)
{
try
{
int id = database::Database::safe_label_to_id(it->first, m_data->labels_minerals);
for (int idc=0; idc<m_data->nb_component; ++idc)
{
m_tot_conc(idc) += m_data->nu_mineral(id, idc)*it->second.first;
m_tot_conc_increase(idc) += m_data->nu_mineral(id, idc)*it->second.second;
}
}
catch(std::invalid_argument& e)
{
int id = database::Database::safe_label_to_id(it->first, m_data->labels_minerals_kinetic);
for (int idc=0; idc<m_data->nb_component; ++idc)
{
m_tot_conc(idc) += m_data->nu_mineral_kinetic(id, idc)*it->second.first;
m_tot_conc_increase(idc) += m_data->nu_mineral_kinetic(id, idc)*it->second.second;
}
}
}
// Do we need to remove some components ?
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);
if (m_model->minerals_to_keep.size() != 0)
{
m_database.minerals_keep_only(m_model->minerals_to_keep);
}
}
micpsolver::MiCPPerformance ReactionPathDriver::one_step(Eigen::VectorXd& x, bool init)
{
ReducedSystemSolver solver;
if (m_current_solution.is_valid())
{
solver = ReducedSystemSolver(m_data, m_tot_conc, m_current_solution);
}
else
{
solver = ReducedSystemSolver(m_data, m_tot_conc);
}
solver.get_options() = m_solver_options;
micpsolver::MiCPPerformance perf = solver.solve(x, init);
if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
{
ERROR << "Failed to solve the system ! (most probably)";
}
m_current_solution = solver.get_solution(x);
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);
}
}
} // end namespace specmicp
c4science · Help