Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/reactmicp/equilibrium_curve/chemistry.cpp
4adc6181a531reactivetransport

Last Change
Raw File

chemistry.cpp
View Options

#include "chemistry.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "utils/log.hpp"
#include <iostream>
namespace specmicp {
namespace reactmicp {
namespace eqcurve {
// EquilibriumCurveSpeciation::
void EquilibriumCurveSpeciation::output()
{
index_t size = m_eqcurve.rows();
if (m_cnt == size) m_eqcurve.conservativeResize(size+500, 4);
m_eqcurve.block(size, 0, m_eqcurve.rows()-size, 4).setZero();
AdimensionalSystemSolutionExtractor sol(current_solution(), m_data, get_units());
m_eqcurve(m_cnt, 0) = sol.total_solid_concentration(m_idc);
m_eqcurve(m_cnt, 1) = sol.total_aqueous_concentration(m_idc);
m_eqcurve(m_cnt, 2) = sol.porosity();
m_eqcurve(m_cnt, 3) = m_eqcurve(m_cnt, 2)>0.92?2.219e-5:1e4*std::exp(9.95*m_eqcurve(m_cnt, 2)-29.08);
std::cout << m_eqcurve(m_cnt, 0) << " \t " << m_eqcurve(m_cnt, 1) << std::endl;
}
Matrix EquilibriumCurveSpeciation::get_equilibrium_curve(scalar_t end_total_concentration, scalar_t delta)
{
m_cnt = 0;
solve_first_problem();
output();
AdimensionalSystemSolutionExtractor sol(current_solution(), m_data, get_units());
index_t nb_steps = std::abs(end_total_concentration - conditions().total_concentrations(m_idc))/(
std::abs(delta)*sol.total_aqueous_concentration(m_idc));
std::cout << nb_steps << std::endl;
m_eqcurve = Matrix::Zero(nb_steps, 4);
m_delta = delta; ///(nb_steps-1);
while (conditions().total_concentrations(m_idc) > end_total_concentration)
{
//std::cout << m_cnt << " - " << conditions().total_concentrations(m_idc) << std::endl;
run_step();
++m_cnt;
}
m_eqcurve.conservativeResize(m_cnt, 4);
return m_eqcurve;
}
void EquilibriumCurveSpeciation::update_problem()
{
Vector diff_conc = get_perturbation();
scalar_t coeff = m_delta;
for (index_t component: m_data->range_component())
{
if (std::abs(conditions().total_concentrations(component)) < 1e-4) {
conditions().total_concentrations(component) = 0.0;
diff_conc(component) = 0.0;
}
if (conditions().total_concentrations(component) != 0 and
std::copysign(1.0, conditions().total_concentrations(component) + coeff*diff_conc(component))
!= std::copysign(1.0, conditions().total_concentrations(component))
)
{
//std::cout << "Component : " << component << std::endl;
coeff = std::copysign(0.9*conditions().total_concentrations(component), m_delta);
coeff /= diff_conc(component);
}
}
//std::cout << "Update : " << std::endl << coeff*diff_conc << std::endl;
conditions().total_concentrations += coeff*diff_conc;
// std::cout << "total concentrations : " << std::endl << conditions().total_concentrations << std::endl;
}
Vector EquilibriumCurveSpeciation::get_perturbation()
{
AdimensionalSystemSolutionExtractor sol(current_solution(), m_data, get_units());
Vector tot_aq_conc(m_data->nb_component);
tot_aq_conc(0) = 0;
for (index_t component: m_data->range_aqueous_component())
{
tot_aq_conc(component) = sol.total_aqueous_concentration(component);
}
//tot_aq_conc *= m_delta/tot_aq_conc.norm();
return tot_aq_conc;
}
void EquilibriumCurveSpeciation::error_handling(std::string msg) const
{
ERROR << msg;
ERROR << " Total concentration " << std::endl << conditions().total_concentrations << std::endl;
throw std::runtime_error(msg);
}
} // end namespace eqcurve
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help