Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
c46a22fd727breactivetransport

Last Change
Raw File

transport_program.cpp
View Options

#include "transport_program.hpp"
#include <iostream>
#define EPS_J 1e-8
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace siasaturated {
void SaturatedDiffusionProgram::number_equations(std::vector<index_t> list_fixed_nodes
)
{
m_ideq = Eigen::VectorXd::Zero(get_ndf()*m_mesh->nnodes());
for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
{
for (index_t component: m_data->range_aqueous_component())
{
m_ideq(get_dof_component(*it, component)) = no_equation;
}
}
index_t neq = 0;
for (index_t dof=0; dof<m_ideq.rows(); ++dof)
{
if (m_ideq(dof) != no_equation)
{
m_ideq(dof) = neq;
++neq;
}
}
m_neq = neq;
}
void SaturatedDiffusionProgram::compute_residuals(
const Vector& displacement,
const Vector& velocity,
Vector& residual
)
{
residual = Eigen::VectorXd::Zero(get_neq());
for (index_t element: m_mesh->range_elements())
{
element_residual_transport(element, displacement, velocity, residual);
}
}
void SaturatedDiffusionProgram::element_residual_transport(index_t element,
const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Eigen::VectorXd& residual)
{
Eigen::VectorXd element_residual(2);
for (index_t component: m_data->range_aqueous_component())
{
element_residual_transport_component(element, component, displacement, velocity, element_residual);
for (index_t en: m_mesh->range_nen())
{
const index_t node = m_mesh->get_node(element, en);
const index_t id = m_ideq(get_dof_component(node, component));
if (id != no_equation) {residual(id) += element_residual(en);}
}
}
}
void SaturatedDiffusionProgram::element_residual_transport_component(index_t element,
index_t component,
const Vector &displacement,
const Vector &velocity,
Vector& element_residual)
{
Eigen::Matrix<scalar_t, 2, 2> mass, jacob;
Eigen::Matrix<scalar_t, 2, 1> evelocity, edisplacement;
scalar_t mass_coeff = -(m_param->density_water() *
m_param->porosity(0) *
m_mesh->crosssection()*m_mesh->get_dx(element)/2);
mass << 1, 0, 0, 1;
mass *= mass_coeff;
scalar_t flux_coeff = -( m_mesh->crosssection() / m_mesh->get_dx(element)
* m_param->porosity(0)
* m_param->density_water()
* m_param->diffusion_coefficient(element)
);
jacob << 1, -1, -1, 1;
jacob *= flux_coeff;
evelocity << velocity(get_dof_component(m_mesh->get_node(element, 0), component)),
velocity(get_dof_component(m_mesh->get_node(element, 1), component));
edisplacement << displacement(get_dof_component(m_mesh->get_node(element, 0), component)),
displacement(get_dof_component(m_mesh->get_node(element, 1), component));
element_residual = mass*evelocity+jacob*edisplacement;
for (index_t en: m_mesh->range_nen())
{
element_residual(en) += (m_mesh->crosssection()*m_mesh->get_dx(element)/2
*external_flux(m_mesh->get_node(element, en), component));
}
}
void SaturatedDiffusionProgram::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
list_triplet_t jacob;
const index_t ncomp = m_data->nb_component-1;
const index_t estimation = m_mesh->nnodes()*(ncomp*m_mesh->nen);
jacob.reserve(estimation);
for (index_t element: m_mesh->range_elements())
{
element_jacobian_transport(element, displacement, velocity, jacob, alphadt);
}
jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
jacobian.setFromTriplets(jacob.begin(), jacob.end());
}
// Transport
// =========
void SaturatedDiffusionProgram::element_jacobian_transport(
index_t element,
Vector& displacement,
Vector& velocity,
list_triplet_t& jacobian,
scalar_t alphadt)
{
for (index_t component: m_data->range_aqueous_component())
{
Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
element_residual_transport_component(element, component, 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 = get_dof_component(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_J*std::abs(tmp_v);
if (h == 0) h = EPS_J;
velocity(dof) = tmp_v + h;
h = velocity(dof) - tmp_v;
displacement(dof) = tmp_d + alphadt*h;
element_residual_transport_component(element, component, 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 = m_ideq(get_dof_component(noder, component));
if (idr == no_equation) continue;
jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
}
}
}
}
void SaturatedDiffusionProgram::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())
{
for (index_t component: m_data->range_aqueous_component())
{
const index_t dof = get_dof_component(node, component);
const index_t id = m_ideq(dof);
if (id == no_equation) continue;
velocity(dof) += lambda*update(id);
}
}
displacement = predictor + alpha_dt*velocity;
}
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help