Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/src/specmicp/kinetics/kinetic_system.cpp
2f65d5a63c08master

Last Change
Raw File

kinetic_system.cpp
View Options

#include "kinetic_system.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include "kinetic_model.hpp"
namespace specmicp {
namespace kinetics {
void KineticSystem::update_total_concentrations(scalar_t dt, const Vector& y)
{
Vector update_minerals = (y - m_variables.moles_minerals())/dt;
RawDatabasePtr database = m_variables.get_database();
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: database->range_component())
{
m_variables.total_concentration(component) +=
database->nu_mineral_kinetic(*it, component)*update_minerals(it);
}
}
}
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_variables.get_database(),
m_variables.total_concentrations(),
m_variables.equilibrium_state(),
options);
variables = m_variables.equilibrium_state().main_variables();
perf = solver.solve(variables);
}
else
{
solver = ReducedSystemSolver(m_variables.get_database(),
m_variables.total_concentrations(),
options);
perf = solver.solve(variables, true);
}
if ((int) perf.return_code < 0)
{
std::cout << "Failed to solve the system ! Err code " << (int) perf.return_code << std::endl;
}
m_variables.update_equilibrium_state(solver.get_solution(variables));
}
} // end namespace kinetics
} // end namespace specmicp
c4science · Help