Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/transport_program.hpp
e5484eae1665yaml/safe

Last Change
Raw File

transport_program.hpp
View Options

#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
#define SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
#include <memory>
#include "database/data_container.hpp"
#include "dfpmsolver/parabolic_program.hpp"
#include "transport_parameters.hpp"
#include "reactmicp/meshes/mesh1d.hpp"
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace siasaturated {
class SaturatedDiffusionProgram: dfpmsolver::ParabolicProgram<SaturatedDiffusionProgram>
{
public:
SaturatedDiffusionProgram(
std::shared_ptr<mesh::Mesh1D> the_mesh,
std::shared_ptr<database::DataContainer> the_data,
std::shared_ptr<SaturatedDiffusionTransportParameters> the_param
):
m_mesh(the_mesh),
m_data(the_data),
m_param(the_param),
m_external_flux(Eigen::VectorXd::Zero(the_mesh->nnodes()*(the_data->nb_component-1)))
{}
const int no_equation = -1;
//! \brief Return the number of equations
int get_neq() const {return m_neq;}
//! \brief Return the number of degrees of freedom per node
int get_ndf() const {return m_data->nb_component-1;}
// number the equations
void number_equations(std::vector<int> list_fixed_nodes);
//! \brief Compute the residuals
void compute_residuals(const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Eigen::VectorXd& residual
);
//! \brief Compute the jacobian
void compute_jacobian(Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity,
Eigen::SparseMatrix<double> &jacobian,
double alphadt);
//! \brief Update the solutions
void update_solution(const Eigen::VectorXd& update,
double lambda,
double alpha_dt,
Eigen::VectorXd& predictor,
Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity);
//! \brief Apply boundary conditions to the velocity vector
//!
//! by default do nothing.
void apply_bc(double dt,
const Eigen::VectorXd& displacement,
Eigen::VectorXd velocity)
{}
//! Return the degree of freedom number of 'component' at 'node'
int get_dof_component(int node, int component) const {
assert(component >=1 and component < m_data->nb_component);
assert(node >= 0 and node < m_mesh->nnodes());
return (component - 1)+get_ndf()*node;
}
//! Return the degree of freedom number of 'component' at 'node'
int get_dof_mineral(int node, int mineral) const {
assert(mineral >= 0 and mineral < m_data->nb_mineral);
assert(node >= 0 and node < m_mesh->nnodes());
return mineral+m_data->nb_mineral*node;
}
//! Return the total mobile concentration
double& get_total_mobile_concentration(int node, int component, Eigen::VectorXd& concentrations) const
{
return concentrations.coeffRef(get_dof_component(node, component));
}
//! Return the total mobile concentration
double get_total_mobile_concentration(int node, int component, const Eigen::VectorXd& concentrations) const
{
return concentrations.coeff(get_dof_component(node, component));
}
//! Return the external flux
Eigen::VectorXd& external_flux() {return m_external_flux;}
//! Return the external flux for 'component' at 'node'
double& external_flux(int node, int component) {return m_external_flux(get_dof_component(node, component));}
private:
void element_residual_transport(
ind_t element,
const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Eigen::VectorXd& residual);
void element_residual_transport_component(
ind_t element,
int component,
const Eigen::VectorXd& displacement,
const Eigen::VectorXd& velocity,
Vector& element_residual);
double nodal_mineral_transient_term_transport(
ind_t node,
ind_t component);
void element_jacobian_transport(
ind_t element,
Eigen::VectorXd& displacement,
Eigen::VectorXd& velocity,
list_triplet_t& jacobian,
double alphadt);
int m_neq;
Eigen::VectorXd m_ideq;
std::shared_ptr<mesh::Mesh1D> m_mesh;
std::shared_ptr<database::DataContainer> m_data;
std::shared_ptr<SaturatedDiffusionTransportParameters> m_param;
Eigen::VectorXd m_external_flux;
};
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
c4science · Help