Page MenuHomec4science
Files F77873456

eqcurve_coupler.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 16, 22:37

eqcurve_coupler.cpp
View Options

#include "eqcurve_coupler.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include "dfpm/1dtransport/diffusion_parameters.hpp"
#include "dfpmsolver/parabolic_driver.hpp"
namespace specmicp {
namespace reactmicp {
namespace eqcurve {
EquilibriumCurveCoupler::EquilibriumCurveCoupler(Matrix& eq_curve,
mesh::Mesh1DPtr the_mesh
, dfpmsolver::ParabolicDriverOptions options):
m_eqcurve(eq_curve),
m_mesh(the_mesh),
m_param(std::make_shared<dfpm::SaturatedDiffusion1DParameters>(m_mesh->nb_nodes())),
m_aqueous_concentrations(the_mesh->nb_nodes()),
m_solid_concentrations(the_mesh->nb_nodes()),
m_transport_program(the_mesh, m_param, {0,}),
m_transport_solver(m_transport_program)
{
m_transport_solver.get_options() = options;
index_t last = m_eqcurve.last();
m_aqueous_concentrations(0) = m_eqcurve.totaq_concentration(last);
m_solid_concentrations(0) = m_eqcurve.totsolid_concentration(last);
m_param->porosity(0) = m_eqcurve.porosity(last);
m_param->diffusion_coefficient(0) = m_eqcurve.diffusion_coefficient(last);
index_t first = m_eqcurve.first();
for (index_t node=1; node<the_mesh->nb_nodes(); ++node)
{
m_aqueous_concentrations(node) = m_eqcurve.totaq_concentration(first);
m_solid_concentrations(node) = m_eqcurve.totsolid_concentration(first);
m_param->porosity(node) = m_eqcurve.porosity(first);
m_param->diffusion_coefficient(node) = m_eqcurve.diffusion_coefficient(first);
}
}
void EquilibriumCurveCoupler::run_step(scalar_t timestep)
{
dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.solve_timestep(timestep, m_aqueous_concentrations);
if (retcode <= dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet)
{
throw std::runtime_error("Cannot solve the transport problem, retcode : "+std::to_string((int) retcode));
}
chemistry_step();
}
void EquilibriumCurveCoupler::chemistry_step()
{
for (index_t node=1; node<m_mesh->nb_nodes(); ++node)
{
index_t index = m_eqcurve.find_point(m_solid_concentrations(node));
//std::cout << index << std::endl;
scalar_t ceq = m_eqcurve.totaq_concentration(index);
//std::cout << m_solid_concentrations(node) << " - diff : " << m_eqcurve.porosity(index)*(ceq-m_aqueous_concentrations(node)) << std::endl;
m_solid_concentrations(node) -= std::min( m_solid_concentrations(node),
m_eqcurve.porosity(index)*(ceq-m_aqueous_concentrations(node)));
m_aqueous_concentrations(node) = ceq;
// ###FIXME -> locate node using previous guess
index = m_eqcurve.find_point(m_solid_concentrations(node));
m_param->porosity(node) = m_eqcurve.porosity(index);
m_param->diffusion_coefficient(node)= m_eqcurve.diffusion_coefficient(index);
}
}
} // end namespace eqcurve
} // end namespace reactmicp
} // end namespace specmicp

Event Timeline

c4science · Help