Homec4science
Diffusion SpecMiCP / ReactMiCP (kinetics) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/transport_neutrality_program.hpp
9a4089e18a62kinetics

Last Change
Raw File

transport_neutrality_program.hpp
View Options

#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSIONNEUTRALITY_TRANSPORTPROGRAM_HPP
#define SPECMICP_REACTMICP_SATURATED_DIFFUSIONNEUTRALITY_TRANSPORTPROGRAM_HPP
#include <memory>
#include "database/data_container.hpp"
#include "dfpmsolver/parabolic_program.hpp"
#include "dfpm/meshes/mesh1dfwd.hpp"
#include "variables.hpp"
#include "transport_neutrality_variables.hpp"
namespace specmicp {
class EquilibriumState;
namespace reactmicp {
namespace systems {
namespace siasaturated {
class SaturatedNeutralityDiffusionTransportParameters;
//! \brief Transport program for saturated diffusion
class SaturatedNeutralityDiffusionProgram: dfpmsolver::ParabolicProgram<SaturatedNeutralityDiffusionProgram>
{
public:
using triplet_t = Eigen::Triplet<scalar_t>; //!< Triplet type for building the transport matrix
using list_triplet_t = std::vector<triplet_t>; //!< List of triplet for building the sparse matrix
SaturatedNeutralityDiffusionProgram(
std::shared_ptr<mesh::Mesh1D> the_mesh,
std::shared_ptr<database::DataContainer> the_data,
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> the_param
);
//! \brief Return the number of equations
index_t get_neq() const {return m_neq;}
//! \brief Return the number of degrees of freedom per node
index_t get_ndf() const {return m_variables.get_ndf();}
index_t get_tot_ndf() const;
// \brief number the equations
void number_equations(std::vector<index_t> list_fixed_nodes);
//! \brief Compute the residuals
void compute_residuals(const Vector &displacement,
const Vector &velocity,
Vector &residual
);
//! \brief Compute the jacobian
void compute_jacobian(Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t> &jacobian,
scalar_t alphadt);
//! \brief Update the solutions
void update_solution(const Vector& update,
scalar_t lambda,
scalar_t alpha_dt,
Vector& predictor,
Vector& displacement,
Vector& velocity);
//! \brief Apply boundary conditions to the velocity vector
//!
//! by default do nothing.
void apply_bc(scalar_t dt,
const Vector& displacement,
Vector velocity)
{}
//! \brief Return the degree of freedom number of 'component' at 'node'
index_t get_dof_component(index_t node, index_t component) const {
return m_variables.get_dof_component(node, component);
}
//! \brief Return the degree of freedom number of 'component' at elemental node 'enode'
index_t get_edof_component(index_t enode, index_t component) const {
return m_variables.get_edof_component(enode, component);
}
index_t id_equation(index_t id_dof) const {
return m_ideq(id_dof);
}
index_t id_equation(index_t node, index_t component) const {
return m_ideq(get_dof_component(node, component));
}
//! \brief Compute total mobile concentrations
void get_total_mobile_concentrations(const Vector& displacement, Vector& total_concentrations);
//! \brief Update concentrations and activity coefficients at node x using an EquilibriumState
void extern_update_nodal_composition(
index_t node,
const EquilibriumState& composition,
Vector& displacement
);
//! \brief Update concentrations and activity coefficients at node x using an EquilibriumState
void extern_update_nodal_composition(
index_t node,
const EquilibriumState& composition
);
//! \brief Return the external flux
Vector& external_flux() {return m_external_flux;}
//! \brief the external flux for 'component' at 'node'
scalar_t& external_flux(index_t node, index_t component) {return m_external_flux(get_dof_component(node, component));}
private:
void copy_displacements_in_variable(
const Vector& displacement
);
void element_residual(
index_t element,
const Vector& velocity,
Vector& eresiduals
);
void element_residual_transport_component_diffusion(
index_t element,
index_t component,
Vector& flux
);
void element_residual_transport_component_nonideal(
index_t element,
index_t component,
Vector& flux,
Vector& potential_factor
);
void element_residual_flux(
index_t element,
Vector& element_residual
);
void element_residual_mass(
index_t element,
const Vector& velocity,
Vector& element_residual
);
void element_residual_mass_fd(
index_t element,
const Vector& velocity,
Vector& element_residual
);
void element_residual_source(
index_t element,
Vector& element_residual
);
void element_jacobian_transport(
index_t element,
Vector& velocity,
list_triplet_t& jacobian,
scalar_t alphadt);
index_t m_neq;
Vector m_ideq;
TransportNeutralityVariables m_variables;
mesh::Mesh1DPtr m_mesh;
database::RawDatabasePtr m_data;
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> m_param;
Vector m_external_flux;
};
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSIONNEUTRALITY_TRANSPORTPROGRAM_HPP
c4science · Help