Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

transport_neutrality_program.hpp
998871ce7287reactivetransport

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

specmicp/src/reactmicp/systems/saturated_diffusion/transport_neutrality_program.hpp998871ce7287reactivetransport

transport_neutrality_program.hppView Options

specmicp/src/reactmicp/systems/saturated_diffusion/transport_neutrality_program.hpp
998871ce7287reactivetransport

transport_neutrality_program.hpp
View Options