Homec4science
Diffusion SpecMiCP / ReactMiCP (configuration) Browse

specmicp/src/reactmicp/systems/saturated_react/transport_program.cpp
3c484d96fa9bconfiguration

Last Change
Raw File

transport_program.cpp
View Options

#include "transport_program.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include "variables.hpp"
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace satdiff {
//class SaturatedDiffusion::
SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
std::vector<index_t> list_fixed_nodes):
m_mesh(variables->get_mesh()),
m_variables(variables),
m_is_in_residual_computation(false)
{
number_equations(list_fixed_nodes);
}
//! \brief Return the number of degrees of freedom per node
inline index_t SaturatedDiffusion::get_ndf() const {return m_variables->nb_component();}
//! \brief Return the total number of degrees of freedom
inline index_t SaturatedDiffusion::get_tot_ndf() const {return m_variables->nb_component()*m_variables->get_mesh()->nb_nodes();}
void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes)
{
m_ideq.resizeLike(m_variables->displacement());
for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
{
for (index_t component=0; component<m_variables->nb_component(); ++component)
{
m_ideq(m_variables->dof_aqueous_concentration(*it, component)) = no_equation;
}
}
index_t neq = 0;
for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
{
for (index_t component=0; component<m_variables->nb_component(); ++component)
{
index_t dof = m_variables->dof_aqueous_concentration(node, component);
if (m_ideq(dof) != no_equation)
{
m_ideq(dof) = neq;
++neq;
}
dof = m_variables->dof_solid_concentration(node, component);
m_ideq(dof) = no_equation;
}
}
m_neq = neq;
}
void SaturatedDiffusion::compute_residuals(const Vector &displacement, const Vector &velocity, Vector &residual)
{
residual = Vector::Zero(get_neq());
m_is_in_residual_computation = true;
for (index_t element: m_mesh->range_elements())
{
for (index_t component=0; component<m_variables->nb_component(); ++component)
{
Eigen::Vector2d element_residual;
element_residual.setZero();
residuals_element_component(element, component, displacement, velocity, element_residual);
for (index_t en=0; en<2; ++en)
{
const index_t node = m_mesh->get_node(element, en);
const index_t id = m_ideq(m_variables->dof_aqueous_concentration(node, component));
if (id != no_equation) {residual(id) += element_residual(en);}
}
}
}
m_is_in_residual_computation = false;
}
void SaturatedDiffusion::residuals_element_component(
index_t element,
index_t component,
const Vector &displacement,
const Vector &velocity,
Eigen::Vector2d &element_residual
)
{
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 index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
const scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) +
0.5/m_variables->diffusion_coefficient(node_1));
scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
* diff_coeff
) ;
const index_t dof_0 = m_variables->dof_aqueous_concentration(node_0, component);
const index_t dof_1 = m_variables->dof_aqueous_concentration(node_1, component);
scalar_t flux = flux_coeff*(displacement(dof_0) - displacement(dof_1));
element_residual(0) = flux;
element_residual(1) = - flux;
if (m_is_in_residual_computation)
{
m_variables->aqueous_concentration(node_0, component,
m_variables->transport_rate()) += element_residual(0);
m_variables->aqueous_concentration(node_1, component,
m_variables->transport_rate()) += element_residual(1);
}
// velocity
element_residual(0) += mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0)
+m_variables->vel_porosity(node_0)*displacement(dof_0));
element_residual(1) += mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1)
+ m_variables->vel_porosity(node_1)*displacement(dof_1));
// external rate
for (index_t en=0; en<2; ++en)
{
element_residual(en) += m_variables->aqueous_concentration(
m_mesh->get_node(element, en), component,
m_variables->chemistry_rate());
}
}
void SaturatedDiffusion::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
dfpm::list_triplet_t jacob;
const index_t ncomp = m_variables->nb_component();
const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
jacob.reserve(estimation);
for (index_t element: m_mesh->range_elements())
{
jacobian_element(element, displacement, velocity, jacob, alphadt);
}
jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
jacobian.setFromTriplets(jacob.begin(), jacob.end());
}
void SaturatedDiffusion::jacobian_element(
index_t element,
Vector& displacement,
Vector& velocity,
dfpm::list_triplet_t& jacobian,
scalar_t alphadt)
{
for (index_t component=0; component<m_variables->nb_component(); ++component)
{
Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero());
residuals_element_component(element, component, displacement, velocity, element_residual_orig);
for (index_t en=0; en<2; ++en)
{
Eigen::Vector2d element_residual(Eigen::Vector2d::Zero());
const index_t node = m_mesh->get_node(element, en);
const index_t dof = m_variables->dof_aqueous_concentration(node, component);
const index_t idc = m_ideq(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-4*eps_jacobian) h = eps_jacobian;
velocity(dof) = tmp_v + h;
h = velocity(dof) - tmp_v;
displacement(dof) = tmp_d + alphadt*h;
residuals_element_component(element, component, displacement, velocity, element_residual);
velocity(dof) = tmp_v;
displacement(dof) = tmp_d;
for (index_t enr=0; enr<2; ++enr)
{
const index_t noder = m_mesh->get_node(element, enr);
const index_t idr = m_ideq(m_variables->dof_aqueous_concentration(noder, component));
if (idr == no_equation) continue;
jacobian.push_back(dfpm::triplet_t(
idr,
idc,
(element_residual(enr) - element_residual_orig(enr))/h
));
}
}
}
}
} // end namespace satdiff
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help