Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85054462
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
Thu, Sep 26, 11:26
Size
10 KB
Mime Type
text/x-c
Expires
Sat, Sep 28, 11:26 (2 d)
Engine
blob
Format
Raw Data
Handle
21123062
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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_ndf(2*variables->nb_component()),
m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes()),
m_mesh(variables->get_mesh()),
m_variables(variables),
m_is_in_residual_computation(false)
{
number_equations(list_fixed_nodes, {}, {0, 1});
}
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_ndf(2*variables->nb_component()),
m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes()),
m_mesh(variables->get_mesh()),
m_variables(variables),
m_is_in_residual_computation(false)
{
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=2; 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=2; 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) // attribute an equation number if it is NOT a slave nor 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=2; 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
Log In to Comment