Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/reduced_system.hpp
1592526f1247reactivetransport

Last Change
Raw File

reduced_system.hpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_REDUCEDSYSTEM_HPP
#define SPECMICP_REDUCEDSYSTEM_HPP
//! \file reduced_system.hpp The MiCP program to solve speciation
#include <memory>
#include "thermodata.hpp"
#include "micpsolver/micpprog.hpp"
//! \namespace specmicp main namespace for the SpecMiCP solver
namespace specmicp {
//! \brief The equilibrium speciation problem
//!
//! Represent the equilibrium speciation problem as a Mixed Complementarity program
//! It is called reduced becaused it does not solve explicitely for the secondary species.
//!
//! Equations
//! ==========
//!
//! Main equations
//! --------------
//!
//! - Water conservation equation : \f[ T_w = m_w \left( \frac{1}{M_w} + \sum_{j=1}^{N_r} \nu_{wj} m_j \right) + \sum_{l=1}^{N_{min}} \nu_{lw} n_l \f]
//! - Component conservation equations : \f[ T_i = m_w \left( m_i + \sum_{j=1}^{N_r} \nu_{ij} m_j \right) + \sum_{l=1}^{N_{min}} \nu_{li} n_l \f]
//! - Mineral equilibrium conditions : \f[ n_l \geq 0, \; 1 - \mathrm{SI}_l \geq 0 , \; \mathrm{and} \; n_l ( 1 - \mathrm{SI}_l ) = 0 \f]
//!
//! Secondary equations
//! -------------------
//!
//! - Secondary aqueous species equilibrium conditions : \f[ m_j = \frac{\prod_{i=2}^{N_c} (\gamma_i m_i)^{\nu_{ji}}}{K_j \gamma_j} \f]
//! - Ionic Strength : \f[ I = \sum_{k=2}^{N_c+N_r} z_k^2 m_k \f]
//! - Activity coefficients : \f[ \log \gamma_k = \frac{-A z_k^2 \sqrt{I}}{1 + B \dot{a}_k \sqrt{I}} + \dot{b}_k I \f]
//!
class ReducedSystem : public micpsolver::MiCPProg<ReducedSystem>
{
public:
const int no_equation = -1; //! The corresponding variables is not a degree of freedom (i.e. is fixed)
//! \brief Create a reduced system
//!
//! \param ptrthermo Shared pointer to the thermodynamic data
//! \param totconc Vector containing the total concentrations (in mol) for each components
//! \param water_eq if true, solve the water conservation equation
ReducedSystem(std::shared_ptr<ThermoData> ptrthermo,
Eigen::VectorXd& totconc,
bool water_eq=true):
m_thermo(ptrthermo),
m_tot_conc(totconc),
m_secondary_conc(ptrthermo->nb_aqueous()),
m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous())
{
number_eq(water_eq);
m_secondary_conc.setZero();
m_loggamma.setZero();
}
// Variables
// ==========
//! \brief Return the total number of variables
int total_variables() const {return m_nb_tot_vars;}
//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
int nb_free_variables() const {return m_nb_free_vars;}
//! \brief Return the number of variables subjected to the complementarity conditions
int nb_complementarity_variables() const {return m_nb_compl_vars;}
// The linear system
// ==================
//! \brief Return the residual for the water conservation equation
double residual_water(const Eigen::VectorXd& x) const;
//! \brief Return the residual for the conservation equation of component (i)
//!
//! For (i==0), water use the method : 'residual_water'
double residual_component(const Eigen::VectorXd& x, int i) const;
//! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$
double residual_mineral(const Eigen::VectorXd& x, int m) const;
//! \brief Return the residuals
void get_residuals(const Eigen::VectorXd& x,
Eigen::VectorXd& residual);
//! \brief Return the jacobian
void get_jacobian(Eigen::VectorXd &x,
Eigen::MatrixXd& jacobian);
//! \brief Return the equation id
int ideq(int i) const {return m_ideq[i];}
//! \brief Return the equation number for conservation of water
int ideq_w() const {return m_ideq[0];}
//! \brief Return the equation number of component 'i'
int ideq_paq(int i) const {
assert(i < m_thermo->nb_component());
return m_ideq[i];
}
//! \brief Return the equation number of mineral 'm'
int ideq_min(int m) const {
assert(m < m_thermo->nb_mineral());
return m_ideq[m+m_thermo->nb_component()];
}
//! \brief Return the mass of water
double mass_water(const Eigen::VectorXd& x) const
{
if (ideq_w() != no_equation) return x(ideq_w());
else return 1.0;
}
//! \brief Compute the ionic strength
void set_ionic_strength(const Eigen::VectorXd& x);
//! \brief Compute the activity coefficients
void compute_log_gamma(const Eigen::VectorXd& x);
//! \brief Compute the activity model, called by the solver at the beginning of an iteration
void hook_start_iteration(const Eigen::VectorXd& x, double norm_residual) {
//compute_log_gamma(x);
not_in_linesearch = true;
if (norm_residual < nb_free_variables())
{
// Use fix point iterations for non-ideality
for (int i=0; i<2; ++i) // 2 seems to reduce the number of iterations, more does not help
{
set_secondary_concentration(x);
compute_log_gamma(x);
}
}
}
void aqueous_tot_conc(const Eigen::VectorXd& x,
Eigen::VectorXd& aq_totconc);
double logIAP_mineral_kinetic(int id_mineral, const Eigen::VectorXd& x);
double logIAP_mineral(int id_mineral, const Eigen::VectorXd& x);
double max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update);
const std::vector<int>& get_non_active_component() {return m_nonactive_component;}
void reasonable_starting_guess(Eigen::VectorXd& x);
void reasonable_restarting_guess(Eigen::VectorXd& x);
private:
// For the normal algorithm
void set_secondary_concentration(const Eigen::VectorXd& x);
// For the tot aq_conc pb
void unsafe_set_secondary_concentration(const Eigen::VectorXd& x);
void number_eq(bool water_equation);
std::shared_ptr<ThermoData> m_thermo;
Eigen::VectorXd m_tot_conc;
int m_nb_tot_vars;
int m_nb_free_vars;
int m_nb_compl_vars;
std::vector<int> m_ideq;
Eigen::VectorXd m_secondary_conc;
double m_ionic_strength;
Eigen::VectorXd m_loggamma;
bool not_in_linesearch;
std::vector<int> m_nonactive_component;
std::vector<bool> m_active_aqueous;
};
} // end namespace specmicp
#endif // SPECMICP_REDUCEDSYSTEM_HPP
c4science · Help