Page MenuHomec4science
Files F82121707

transport_program.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Sep 9, 17:15

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),
m_ndf(2*variables->nb_component()),
m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes())
{
number_equations(list_fixed_nodes, {}, {0,});
}
SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
std::vector<index_t> list_fixed_nodes,
std::map<index_t, index_t> list_slave_nodes,
std::vector<index_t> list_immobile_components):
m_mesh(variables->get_mesh()),
m_variables(variables),
m_is_in_residual_computation(false),
m_ndf(2*variables->nb_component()),
m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes())
{
number_equations(list_fixed_nodes, list_slave_nodes, list_immobile_components);
}
void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes,
std::map<index_t, index_t> list_slave_nodes,
std::vector<index_t> list_immobile_components)
{
m_ideq.resizeLike(m_variables->displacement());
m_ideq.setZero();
// flag fixed nodes
for (index_t node: list_fixed_nodes)
{
for (index_t component=1; component<m_variables->nb_component(); ++component)
{
m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
}
}
// flag slaves nodes
// we flag them by making their ideq more negative than no_equation
for (auto slave_pair: list_slave_nodes)
{
for (index_t component=1; component<m_variables->nb_component(); ++component)
{
m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) = no_equation-1;
}
}
// set equation numbers
index_t neq = 0;
for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
{
for (index_t component: list_immobile_components)
{
m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
}
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) // don't attribute an equation number if it is a slave or a fixed node
{
m_ideq(dof) = neq;
++neq;
}
dof = m_variables->dof_solid_concentration(node, component);
m_ideq(dof) = no_equation;
}
}
// slave nodes
// attribute the correct equation number
for (auto slave_pair: list_slave_nodes)
{
for (index_t component=1; component<m_variables->nb_component(); ++component)
{
m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) =
m_ideq(m_variables->dof_aqueous_concentration(slave_pair.second, component));
}
}
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=1; 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);
// diffusion
scalar_t flux_diffusion = flux_coeff*(displacement(dof_0) - displacement(dof_1));
element_residual(0) = flux_diffusion;
element_residual(1) = - flux_diffusion;
// advection
if (m_variables->fluid_velocity(element) != 0.0)
{
scalar_t flux_advection = (m_mesh->get_face_area(element))
*m_variables->fluid_velocity(element);
if (m_variables->fluid_velocity(element) > 0)
{
flux_advection *= (displacement(dof_0) - displacement(dof_1));
element_residual(1) += flux_advection;
}
else
{
flux_advection *= (displacement(dof_1) - displacement(dof_0));
element_residual(0) -= flux_advection;
}
}
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
element_residual(0) += mass_coeff_0*m_variables->solid_concentration(node_0, component,
m_variables->chemistry_rate());
element_residual(1) += mass_coeff_1*m_variables->solid_concentration(node_1, 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=1; 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
));
}
}
}
}
//! \brief Update the solutions
void SaturatedDiffusion::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=1; component< m_variables->nb_component(); ++component)
{
const index_t dof = m_variables->dof_aqueous_concentration(node, component);
const index_t id = m_ideq(dof);
if (id == no_equation) continue;
velocity(dof) += lambda*update(id);
}
}
//displacement = m_variables->predictor() + alpha_dt*velocity;
displacement = predictor + alpha_dt*velocity;
}
} // end namespace satdiff
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp

Event Timeline

c4science · Help