Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86668420
chemistry.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Oct 7, 22:09
Size
3 KB
Mime Type
text/x-c
Expires
Wed, Oct 9, 22:09 (2 d)
Engine
blob
Format
Raw Data
Handle
21464821
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Event Timeline
Log In to Comment