Page MenuHomec4science
Files F76649220

diffusion.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 9, 14:14

diffusion.cpp
View Options

#include "diffusion.hpp"
#include "diffusion_parameters.hpp"
#include "dfpm/meshes/mesh1d.hpp"
namespace specmicp {
namespace dfpm {
// SaturatedDiffusion1D::
SaturatedDiffusion1D::SaturatedDiffusion1D(
mesh::Mesh1DPtr the_mesh,
std::shared_ptr<SaturatedDiffusion1DParameters> parameters,
std::vector<index_t> list_bcs
):
m_tot_ndf(the_mesh->nb_nodes()),
m_mesh(the_mesh),
m_param(parameters),
m_internal_flow(Vector::Zero(m_tot_ndf)),
m_external_flow(Vector::Zero(m_tot_ndf))
{
number_equations(list_bcs);
}
void SaturatedDiffusion1D::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 SaturatedDiffusion1D::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())
{
if (id_equation(m_mesh->get_node(element, enode)) == no_equation) continue;
residual(id_equation(m_mesh->get_node(element, enode))) += elem_residuals(enode);
}
}
}
void SaturatedDiffusion1D::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, edisplacement;
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 porosity = (m_param->porosity(node_0) +
m_param->porosity(node_1) )/2.0;
const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_coefficient(node_0) +
0.5/m_param->diffusion_coefficient(node_1));
const index_t dof_0 = node_0;
const index_t dof_1 = node_1;
scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
// * porosity
* diff_coeff
);
jacob << 1.0, -1.0, -1.0, 1.0;
jacob *= flux_coeff;
evelocity << velocity(dof_0)*
mass_coeff*m_param->porosity(node_0),
velocity(dof_1)*
mass_coeff*m_param->porosity(node_1);
edisplacement << displacement(dof_0),
displacement(dof_1);
elem_residuals = jacob*edisplacement;
m_internal_flow(dof_0) += elem_residuals(0);
m_internal_flow(dof_1) += elem_residuals(1);
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 SaturatedDiffusion1D::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
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());
}
void SaturatedDiffusion1D::element_jacobian(
index_t element,
Vector& displacement,
Vector& velocity,
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-4*eps_jacobian) 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(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
}
}
}
void SaturatedDiffusion1D::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 dfpm
} // end namespace specmicp

Event Timeline

c4science · Help