Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/variables.hpp
e5484eae1665reactmicp/unsaturated

Last Change
Raw File

variables.hpp
View Options

#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP
#define SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP
#include "reactmicp/systems/secondary_variables/secondary.hpp"
#include "specmicp/equilibrium_data.hpp"
#include "reactmicp/meshes/mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace siasaturated {
//! Variables for the SIASaturatedReactiveTransportSolver
class SIASaturatedVariables: public SecondaryVariables<SIASaturatedVariables>
{
public:
SIASaturatedVariables(
std::shared_ptr<specmicp::database::DataContainer> thedatabase,
std::shared_ptr<mesh::Mesh1D> themesh,
std::shared_ptr<SaturatedDiffusionTransportParameters> theparam
):
SecondaryVariables<SIASaturatedVariables>(themesh->nnodes(), thedatabase),
m_data(thedatabase),
m_mesh(themesh),
m_param(theparam),
m_total_mobile_concentration(
Eigen::VectorXd((thedatabase->nb_component-1)* themesh->nnodes())),
m_speciation_variables(
Eigen::MatrixXd(thedatabase->nb_component+thedatabase->nb_mineral, themesh->nnodes()))
{}
int get_dof_total_concentration(int node, int component) const {
assert(component >= 1 and component < m_data->nb_component);
assert(node >= 0 and node < m_speciation_variables.cols());
return (component+offset_total_component())+(m_data->nb_component+offset_total_component())*node;
}
//! Return the total mobile concentration of 'component' at 'node'
double total_mobile_concentration(int node, int component) const {
return m_total_mobile_concentration.coeff(get_dof_total_concentration(node, component));
}
//! Return a reference to the total mobile concentration of 'component' at 'node'
double& total_mobile_concentration(int node, int component) {
return m_total_mobile_concentration.coeffRef(get_dof_total_concentration(node, component));
}
Eigen::VectorXd& total_mobile_concentrations() {
return m_total_mobile_concentration;
}
Eigen::Block<Eigen::VectorXd, Eigen::Dynamic, 1> total_mobile_concentrations(int node) {
return m_total_mobile_concentration.segment(node*(m_data->nb_component+offset_total_component()),
(m_data->nb_component+offset_total_component()));
}
//! Return the concentration of 'component' at 'node'
double component_concentration(int node, int component) const {
assert(component >= 1 and component < m_data->nb_component);
assert(node >= 0 and node < m_speciation_variables.cols());
return pow10(m_speciation_variables.coeff(component+offset_component(), node));
}
//! Return the concentration of 'component' at 'node'
double log_component_concentration(int node, int component) const {
assert(component >= 1 and component < m_data->nb_component);
assert(node >= 0 and node < m_speciation_variables.cols());
return m_speciation_variables.coeff(component+offset_component(), node);
}
//! Return a reference to the concentration of 'component' at 'node'
double& log_component_concentration(int node, int component) {
assert(component >= 1 and component < m_data->nb_component);
assert(node >= 0 and node < m_speciation_variables.cols());
return m_speciation_variables.coeffRef(component+offset_component(), node);
}
//! Return the number of moles of 'mineral' at 'node'
double mineral_amount(int node, int mineral) const {
assert(mineral >= 0 and mineral < m_data->nb_mineral);
assert(node >= 0 and node < m_speciation_variables.cols());
return m_speciation_variables.coeff(offset_mineral()+mineral, node);
}
//! Return a reference to the number of moles of 'mineral' at 'node'
double& mineral_amount(int node, int mineral) {
assert(mineral >= 0 and mineral < m_data->nb_mineral);
assert(node >= 0 and node < m_speciation_variables.cols());
return m_speciation_variables.coeffRef(offset_mineral()+mineral, node);
}
//! Reference to the speciation variables
Eigen::MatrixXd& speciation_variables() {
return m_speciation_variables;
}
//! Eigen Expression to the speciation variables
Eigen::MatrixXd::ColXpr speciation_variables(int node) {
return m_speciation_variables.col(node);
}
//! Return the chemical composition of a node - needed to initialize the speciation solver
EquilibriumState equilibrium_composition(int node) {
return EquilibriumState(m_speciation_variables.col(node),
m_secondary_concentration.col(node),
m_secondary_variables.col(node),
m_ionic_strength(node),
m_data
);
}
//! Update the composition at 'node' using the EquilibriumState instance 'composition'
void update_composition(int node, const EquilibriumState &composition);
//! Return the total moles number for 'component' at 'node' from the speciation
double immobile_total_amount(int node, int component);
//! Return the mass of water at 'node'
double mass_water(int node);
//! Compute the total moles number at 'node' using the transport mobile concentration
void nodal_update_total_amount(int node, Eigen::VectorXd& total_amount);
// Used in the update FEM->speciation
double nodal_component_update_total_amount(int node, int component);
private:
// where the component informations are stored
int offset_total_component() const {return -1;}
// where the component informations are stored
int offset_component() const {return 0;}
// where the mineral informations are stored
int offset_mineral() const {return m_data->nb_component+offset_component();}
std::shared_ptr<specmicp::database::DataContainer> m_data;
std::shared_ptr<mesh::Mesh1D> m_mesh;
std::shared_ptr<SaturatedDiffusionTransportParameters> m_param;
Eigen::VectorXd m_total_mobile_concentration;
Eigen::MatrixXd m_speciation_variables;
};
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP
c4science · Help