kinetic_system_euler_solver.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

	#include "kinetic_system_euler_solver.hpp"

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

	#include <functional>

	namespace specmicp {
	namespace kinetics {


	void AdimKineticSystemEulerSolver::solve(scalar_t dt, scalar_t total)
	{
	double t = 0.0;
	Vector y = m_system.variables().concentration_minerals();
	Vector dydx(y.rows());
	dydx.setZero();

	while (t < total)
	{
	m_system.rhs(t, y, dydx, get_options().speciation_options);
	y += dt*dydx;

	m_system.update_to_new_initial_condition(y, dt);
	t += dt;
	if (t+dt > total) dt =(total - t);
	}
	m_system.update_total_concentrations(y);
	m_system.compute_equilibrium(m_system.constraints(), get_options().speciation_options);

	}

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