Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/src/specmicp/kinetics/kinetic_system.cpp
8ec8f50ef0fayaml/safe

Last Change
Raw File

kinetic_system.cpp
View Options

#include "kinetic_system.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include "kinetic_model.hpp"
#include "utils/log.hpp"
namespace specmicp {
namespace kinetics {
void KineticSystem::update_total_concentrations(const Vector& y)
{
//Vector y_temp = y.cwiseMax(Vector::Zero(y.rows()));
Vector update_minerals = (y - m_variables.moles_minerals());
m_variables.total_concentrations() = m_variables.total_concentrations_initial();
for (auto it=m_model->species_begin(); it!=m_model->species_end(); ++it)
{
for (index_t component: m_data->range_component())
{
if (m_data->nu_mineral_kinetic(*it, component) == 0.0) continue;
m_variables.total_concentration(component) -=
m_data->nu_mineral_kinetic(*it, component)*update_minerals(it-m_model->species_begin());
}
}
}
void KineticSystem::update_to_new_initial_condition(const Vector& y, scalar_t dt)
{
update_total_concentrations(y);
m_variables.rate_components() += (
m_variables.total_concentrations()
- m_variables.total_concentrations_initial()
)/dt;
m_variables.moles_minerals() = y;
m_variables.total_concentrations_initial() = m_variables.total_concentrations();
}
void KineticSystem::compute_rates(scalar_t t, const Vector& y, Vector& dydt)
{
m_model->compute_rate(t, y, m_variables, dydt);
}
void KineticSystem::compute_equilibrium(ReducedSystemSolverOptions options)
{
ReducedSystemSolver solver;
Vector variables;
specmicp::micpsolver::MiCPPerformance perf;
if (m_variables.equilibrium_state().is_valid())
{
solver = ReducedSystemSolver(m_data,
m_variables.total_concentrations(),
m_variables.equilibrium_state(),
options);
variables = m_variables.equilibrium_state().main_variables();
perf = solver.solve(variables);
}
else
{
solver = ReducedSystemSolver(m_data,
m_variables.total_concentrations(),
options);
perf = solver.solve(variables, true);
}
//std::cout << variables << std::endl;
if ((int) perf.return_code < 0)
{
ERROR << "Failed to solve the system ! Error code " << (int) perf.return_code;
}
m_variables.update_equilibrium_state(solver.get_solution(variables));
}
} // end namespace kinetics
} // end namespace specmicp
c4science · Help