Homec4science
Diffusion SpecMiCP / ReactMiCP (single) Browse

specmicp/src/reactmicp/systems/secondary_variables/secondary.inl
c46a22fd727bsingle

Last Change
Raw File

secondary.inl
View Options

#include "secondary.hpp"
#include "physics/laws.hpp"
namespace specmicp {
namespace reactmicp {
// ================================== //
// //
// Secondary variables wibbly-timey //
// //
// ================================== //
//! \brief Set secondary concentrations
template <class Derived>
void SecondaryVariables<Derived>::compute_secondary_concentrations()
{
for (index_t node=0; node<m_secondary_concentration.cols(); ++node)
{
nodal_secondary_concentrations(node);
}
}
//! \brief Set secondary concentrations for node 'node'
template <class Derived>
void SecondaryVariables<Derived>::nodal_secondary_concentrations(index_t node)
{
for (int species: m_data->range_aqueous())
{
scalar_t sum = -m_data->logk_aqueous(species) - log_gamma_secondary(node, species);
for (int component: m_data->range_aqueous_component())
{
sum += m_data->nu_aqueous(species, component)*(
log_gamma_component(node, component)
+ log_component_concentration(node, component)
);
}
secondary_concentration(node, species) = pow10(sum);
}
}
template <class Derived>
double SecondaryVariables<Derived>::norm_secondary_variables()
{
return m_secondary_variables.block(
1, 0,
m_data->nb_aqueous+m_data->nb_component,
m_secondary_concentration.cols()
).norm();
}
template <class Derived>
double SecondaryVariables<Derived>::nodal_norm_secondary_variables(index_t node)
{
return m_secondary_variables.block(
1, node,
m_data->nb_aqueous+m_data->nb_component,
1
).norm();
}
//! \brief Solve for secondary variables
template <class Derived>
int SecondaryVariables<Derived>::solve_secondary_variables()
{
bool may_have_converged = false;
scalar_t previous_norm = norm_secondary_variables();
for (int i=0; i<6; ++i)
{
compute_secondary_concentrations();
compute_secondary_variables();
double new_norm = norm_secondary_variables();
if (std::abs(previous_norm - new_norm)/previous_norm < 1e-6) {may_have_converged=true; break;}
previous_norm = new_norm;
}
return may_have_converged;
}
//! \brief Solve for secondary variables
template <class Derived>
int SecondaryVariables<Derived>::nodal_solve_secondary_variables(index_t node)
{
bool may_have_converged = false;
scalar_t previous_norm = nodal_norm_secondary_variables(node);
for (int i=0; i<6; ++i)
{
nodal_secondary_concentrations(node);
nodal_secondary_variables(node);
double new_norm = nodal_norm_secondary_variables(node);
if (std::abs(previous_norm - new_norm)/previous_norm < 1e-6) {may_have_converged=true; break;}
previous_norm = new_norm;
}
return may_have_converged;
}
//! \brief Compute secondary variables
template <class Derived>
void SecondaryVariables<Derived>::compute_secondary_variables()
{
for (index_t node=0; node<m_secondary_concentration.cols(); ++node)
{
nodal_secondary_variables(node);
}
}
template <class Derived>
void SecondaryVariables<Derived>::nodal_secondary_variables(index_t node)
{
nodal_ionic_strength(node);
nodal_log_gamma(node);
}
template <class Derived>
void SecondaryVariables<Derived>::compute_ionic_strength()
{
for (index_t node=0; node<m_secondary_concentration.cols(); ++node)
{
nodal_ionic_strength(node);
}
}
//! \brief compute the ionic strength for node 'node'
template <class Derived>
void SecondaryVariables<Derived>::nodal_ionic_strength(index_t node)
{
scalar_t ionics = 0;
for (int component: m_data->range_aqueous_component())
{
if (m_data->param_aq(component, 0) == 0) continue;
ionics += std::pow(m_data->param_aq(component, 0), 2)*component_concentration(node, component);
}
for (int aqueous: m_data->range_aqueous())
{
index_t idaq = m_data->nb_component + aqueous;
if (m_data->param_aq(idaq, 0) == 0) continue;
ionics += std::pow(m_data->param_aq(idaq, 0), 2)*secondary_concentration(node, aqueous);
}
ionic_strength(node) = ionics/2;
}
template <class Derived>
void SecondaryVariables<Derived>::compute_log_gamma()
{
for (index_t node=0; node<m_secondary_concentration.cols(); ++node)
{
nodal_log_gamma(node);
}
}
//! \brief compute the logarithm of the activity coefficients for node 'node'
template <class Derived>
void SecondaryVariables<Derived>::nodal_log_gamma(index_t node)
{
nodal_log_gamma_component(node);
nodal_log_gamma_aqueous(node);
}
//! \brief compute the logarithm of the activity coefficients for node 'node'
template <class Derived>
void SecondaryVariables<Derived>::nodal_log_gamma_component(index_t node)
{
const scalar_t sqrti = std::sqrt(ionic_strength(node));
for (int component: m_data->range_aqueous_component())
{
log_gamma_component(node, component) = laws::extended_debye_huckel(ionic_strength(node), sqrti,
m_data->param_aq(component, 0),
m_data->param_aq(component, 1),
m_data->param_aq(component, 2));
}
}
//! \brief compute the logarithm of the activity coefficients for node 'node'
template <class Derived>
void SecondaryVariables<Derived>::nodal_log_gamma_aqueous(index_t node)
{
const scalar_t sqrti = std::sqrt(ionic_strength(node));
for (int aqueous: m_data->range_aqueous())
{
index_t idaq = m_data->nb_component+aqueous;
log_gamma_secondary(node, aqueous) = laws::extended_debye_huckel(ionic_strength(node), sqrti,
m_data->param_aq(idaq, 0),
m_data->param_aq(idaq, 1),
m_data->param_aq(idaq, 2));
}
}
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help