kinetic_system_solver.cpp
522d1d188e48specmicp_bin

kinetic_system_solver.cpp
View Options

			#include "kinetic_system_solver.hpp"

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

			#include <functional>

			namespace specmicp {
			namespace kinetics {


			void AdimKineticSystemSolver::solve(scalar_t dt, scalar_t total)
			{
			double t = 0.0;
			Vector y = m_system.variables().concentration_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(&AdimKineticSystem::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::DormandPrinceStep<Eigen::Dynamic> driver = odeint::get_dormand_prince_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(m_system.constraints(), get_options().speciation_options);

			}

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