Page MenuHomec4science
Files F92302067

eqcurve_solid_transport.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 19, 05:42

eqcurve_solid_transport.cpp
View Options

#include "eqcurve_solid_transport.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include <iostream>
namespace specmicp {
namespace reactmicp {
namespace eqcurve {
// SolidDiffusion::
SolidDiffusion::SolidDiffusion(
mesh::Mesh1DPtr the_mesh,
const Matrix& eq_curve,
std::vector<index_t> list_bcs
):
m_tot_ndf(the_mesh->nb_nodes()),
m_mesh(the_mesh),
m_eqcurve(eq_curve),
m_internal_flow(Vector::Zero(m_tot_ndf)),
m_external_flow(Vector::Zero(m_tot_ndf)),
m_in_jac(false)
{
number_equations(list_bcs);
}
void SolidDiffusion::number_equations(std::vector<index_t> list_bcs)
{
m_id_equations = Eigen::VectorXi::Zero(m_tot_ndf);
for (auto it=list_bcs.begin(); it!=list_bcs.end(); ++it)
{
m_id_equations(*it) = no_equation;
}
index_t neq = 0;
for (index_t node: m_mesh->range_nodes())
{
if (m_id_equations(node) == no_equation) continue;
m_id_equations(node) = neq;
++ neq;
}
m_neq = neq;
}
void SolidDiffusion::compute_residuals(
const Vector& displacement,
const Vector& velocity,
Vector& residual
)
{
m_internal_flow.setZero();
residual.resize(get_neq());
residual.setZero();
for (index_t element: m_mesh->range_elements())
{
Vector elem_residuals(2);
elem_residuals.setZero();
element_residuals(element, displacement, velocity, elem_residuals);
for (index_t enode: m_mesh->range_nen())
{
const index_t id_eq = id_equation(m_mesh->get_node(element, enode));
if (id_eq == no_equation) continue;
residual(id_eq) += elem_residuals(enode);
}
}
//std::cout << "residual : " << std::endl << residual << std::endl;
//std::cout << "flow : " << std::endl << m_internal_flow << std::endl;
}
void SolidDiffusion::element_residuals(
index_t element,
const Vector& displacement,
const Vector& velocity,
Vector& elem_residuals
)
{
Eigen::Matrix<scalar_t, 2, 2> jacob;
Eigen::Matrix<scalar_t, 2, 1> evelocity, econc;
//scalar_t mass_coeff = -(m_mesh->get_volume_element(element)/2.0);
const index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0);
const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1);
const scalar_t sc_0 = displacement(node_0);
const scalar_t sc_1 = displacement(node_1);
const index_t index_0 = m_eqcurve.find_point(sc_0);
const index_t index_1 = m_eqcurve.find_point(sc_1);
//std::cout << element << " # " << index_0 << " - " << index_1 << std::endl;
const scalar_t cc_0 = m_eqcurve.interpolate(index_0, sc_0, 1);
//std::max(m_eqcurve.totaq_concentration(m_eqcurve.last()), m_eqcurve.interpolate(index_0, sc_0, 1));
const scalar_t cc_1 = m_eqcurve.interpolate(index_1, sc_1, 1);
//std::max(m_eqcurve.totaq_concentration(m_eqcurve.last()), m_eqcurve.interpolate(index_1, sc_1, 1));
const scalar_t porosity_0 = m_eqcurve.interpolate(index_0, sc_0, 2);
const scalar_t porosity_1 = m_eqcurve.interpolate(index_1, sc_1, 2);
//const scalar_t diffcoeff_0 = m_eqcurve.interpolate(index_0, sc_0, 3);
//const scalar_t diffcoeff_1 = m_eqcurve.interpolate(index_1, sc_1, 3);
const scalar_t porosity = ( porosity_0
+ porosity_1
)/2.0;
const scalar_t diff_coeff = porosity>0.92?2.219e-5:1e4*std::exp(9.95*porosity-29.08);
// const scalar_t diff_coeff = 1.0/(0.5/diffcoeff_0 +
// 0.5/diffcoeff_1);
scalar_t flux_coeff = ( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
//* porosity
* diff_coeff
);
// if (m_eqcurve.slope(index_1, 1) != 0)
// {
// std::cout << element << " # " << m_eqcurve.slope(index_0, 1) << " - " << m_eqcurve.slope(index_0, 2)
// << " # " << m_eqcurve.slope(index_1, 1) << " - " << m_eqcurve.slope(index_1, 2) << std::endl;
// }
evelocity << (+ m_eqcurve.slope(index_0, 1)*porosity_0
+ m_eqcurve.slope(index_0, 2)*cc_0
+ 1.0)*mass_coeff_0*velocity(node_0),
(+ m_eqcurve.slope(index_1, 1)*porosity_1
+ m_eqcurve.slope(index_1, 2)*cc_1
+ 1.0)*mass_coeff_1*velocity(node_1);
jacob << 1.0, -1.0, -1.0, 1.0;
jacob *= flux_coeff;
econc << cc_0,
cc_1;
// if (element == 0)
// std::cout << econc << std::endl;
elem_residuals += jacob*econc;
m_internal_flow(node_0) += elem_residuals(0);
m_internal_flow(node_1) += elem_residuals(1);
elem_residuals += evelocity;
// for (index_t en: m_mesh->range_nen())
// {
// elem_residuals(en) += evelocity(en);
// elem_residuals(en) += (m_mesh->get_volume_element(element)/2.0
// *external_flow(m_mesh->get_node(element, en)));
// }
}
void SolidDiffusion::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
m_in_jac = true;
dfpm::list_triplet_t jacob;
const index_t estimation = m_mesh->nb_nodes()*(m_mesh->nen);
jacob.reserve(estimation);
for (index_t element: m_mesh->range_elements())
{
element_jacobian(element, displacement, velocity, jacob, alphadt);
}
jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
jacobian.setFromTriplets(jacob.begin(), jacob.end());
m_in_jac = false;
}
void SolidDiffusion::element_jacobian(
index_t element,
Vector& displacement,
Vector& velocity,
dfpm::list_triplet_t& jacobian,
scalar_t alphadt)
{
Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
element_residuals(element, displacement, velocity, element_residual_orig);
for (index_t en: m_mesh->range_nen())
{
Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
const index_t node = m_mesh->get_node(element, en);
const index_t dof = node;
const index_t idc = id_equation(dof);
if (idc == no_equation) continue;
const scalar_t tmp_v = velocity(dof);
const scalar_t tmp_d = displacement(dof);
scalar_t h = eps_jacobian*std::abs(tmp_v);
if (h < 1e-6) h = eps_jacobian;
velocity(dof) = tmp_v + h;
h = velocity(dof) - tmp_v;
displacement(dof) = tmp_d + alphadt*h;
element_residuals(element, displacement, velocity, element_residual);
velocity(dof) = tmp_v;
displacement(dof) = tmp_d;
for (index_t enr: m_mesh->range_nen())
{
const index_t noder = m_mesh->get_node(element, enr);
const index_t idr = id_equation(noder);
if (idr == no_equation) continue;
jacobian.push_back(dfpm::triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
}
}
}
void SolidDiffusion::update_solution(const Vector &update,
scalar_t lambda,
scalar_t alpha_dt,
Vector &predictor,
Vector &displacement,
Vector &velocity
)
{
for (index_t node: m_mesh->range_nodes())
{
const index_t id = id_equation(node);
if (id == no_equation) continue;
velocity(node) += lambda*update(id);
}
displacement = predictor + alpha_dt*velocity;
}
} // end namespace eqcurve
} // end namespace reactmicp
} // end namespace specmicp

Event Timeline

c4science · Help