kinetic_system_solver.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

	#include "kinetic_system_solver.hpp"

	#include "kinetic_model.hpp"
	#include "odeint/runge_kutta_step.hpp"

	#include <functional>

	namespace specmicp {
	namespace kinetics {


	void KineticSystemSolver::solve(scalar_t dt, scalar_t total)
	{
	double t = 0.0;
	Vector y = m_system.variables().moles_minerals();
	Vector dydx(y.rows());
	dydx.setZero();
	m_system.rhs(t, y, dydx, get_options().speciation_options);


	odeint::rhs_f<Eigen::Dynamic> rhs = std::bind(std::mem_fn(&KineticSystem::rhs),
	m_system,
	std::placeholders::_1,
	std::placeholders::_2,
	std::placeholders::_3,
	get_options().speciation_options);
	odeint::CashKarpStep<Eigen::Dynamic> driver = odeint::get_cash_karp_step<Eigen::Dynamic>(rhs);

	odeint::StepLength stepl(dt, total);
	int cnt = 0;
	while (t < total)
	{
	driver.rk_step(y, dydx, t, stepl);
	m_system.update_to_new_initial_condition(y, stepl.did);
	cnt +=1;
	scalar_t remaining = stepl.remaining(t);
	if (remaining == 0.0) break;
	if (stepl.next > remaining) stepl.next = remaining;
	stepl.get_next();
	m_current_dt = stepl.test;

	}
	m_system.update_total_concentrations(y);
	m_system.compute_equilibrium(get_options().speciation_options);

	}

	} // end namespace kinetics
	} // end namespace specmicp