Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
83e1af765087reactmicp/unsaturated

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<int> 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 (int component: m_data->range_aqueous_component())
{
m_ideq(get_dof_component(*it, component)) = no_equation;
}
}
int neq = 0;
for (int 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 Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Eigen::VectorXd& residual
)
{
residual = Eigen::VectorXd::Zero(get_neq());
for (ind_t element: m_mesh->range_elements())
{
element_residual_transport(element, displacement, velocity, residual);
}
}
void SaturatedDiffusionProgram::element_residual_transport(
ind_t element,
const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Eigen::VectorXd& residual)
{
Eigen::VectorXd element_residual(2);
for (ind_t component: m_data->range_aqueous_component())
{
element_residual_transport_component(element, component, displacement, velocity, element_residual);
for (int en: m_mesh->range_nen())
{
const ind_t node = m_mesh->get_node(element, en);
const ind_t id = m_ideq(get_dof_component(node, component));
if (id != no_equation) {residual(id) += element_residual(en);}
}
}
}
void SaturatedDiffusionProgram::element_residual_transport_component(
ind_t element,
int component,
const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Vector& element_residual)
{
Eigen::Matrix2d mass, jacob;
Eigen::Vector2d evelocity, edisplacement;
double 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;
double 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;
//std::cout << "elem " << std::endl << element_residual << std::endl;
for (int 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));
}
//std::cout << "elem " << std::endl << element_residual << std::endl;
}
void SaturatedDiffusionProgram::compute_jacobian(
Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity,
Eigen::SparseMatrix<double>& jacobian,
double alphadt
)
{
list_triplet_t jacob;
const ind_t ncomp = m_data->nb_component-1;
const ind_t estimation = m_mesh->nnodes()*(ncomp*m_mesh->nen);
jacob.reserve(estimation);
for (ind_t element: m_mesh->range_elements())
{
element_jacobian_transport(element, displacement, velocity, jacob, alphadt);
}
jacobian = SparseMatrix(get_neq(), get_neq());
jacobian.setFromTriplets(jacob.begin(), jacob.end());
}
// Transport
// =========
void SaturatedDiffusionProgram::element_jacobian_transport(
ind_t element,
Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity,
list_triplet_t& jacobian,
double alphadt)
{
for (int 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 (int en: m_mesh->range_nen())
{
Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
const ind_t node = m_mesh->get_node(element, en);
const int dof = get_dof_component(node, component);
const int idc = m_ideq(dof);
if (idc == no_equation) continue;
const double tmp_v = velocity(dof);
const double tmp_d = displacement(dof);
double 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 (int enr: m_mesh->range_nen())
{
const ind_t noder = m_mesh->get_node(element, enr);
const ind_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 Eigen::VectorXd& update,
double lambda,
double alpha_dt,
Eigen::VectorXd& predictor,
Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity
)
{
for (int node: m_mesh->range_nodes())
{
for (int component: m_data->range_aqueous_component())
{
const int dof = get_dof_component(node, component);
const int id = m_ideq(dof);
if (id == no_equation) continue;
velocity(dof) += update(id);
}
}
displacement = predictor + alpha_dt*velocity;
}
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help