Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71804117
transport_program.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jul 13, 05:26
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Jul 15, 05:26 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
18996293
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Event Timeline
Log In to Comment