Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/reduced_system.hpp
ce5e17f2cd94reactivetransport

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 specmicp.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 {
class ReducedSystem : public micpsolver::MiCPProg<ReducedSystem>
{
public:
const int no_equation = -1;
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
int total_variables() const {return m_nb_tot_vars;}
int nb_free_variables() const {return m_nb_free_vars;}
int nb_complementarity_variables() const {return m_nb_compl_vars;}
// Residual
double residual_water(const Eigen::VectorXd& x) const;
double residual_component(const Eigen::VectorXd& x, int i) const;
double residual_equilibrium_aqueous(const Eigen::VectorXd& x, int j) const;
double residual_mineral(const Eigen::VectorXd& x, int m) const;
void get_residuals(const Eigen::VectorXd& x,
Eigen::VectorXd& residual);
void get_jacobian(Eigen::VectorXd &x,
Eigen::MatrixXd& jacobian);
void init_secondary_conc(Eigen::VectorXd& x);
int ideq(int i) const {return m_ideq[i];}
int ideq_w() const {return m_ideq[0];}
int ideq_paq(int i) const {
assert(i < m_thermo->nb_component());
return m_ideq[i];
}
int ideq_min(int m) const {
assert(m < m_thermo->nb_mineral());
return m_ideq[m+m_thermo->nb_component()];
}
double mass_water(const Eigen::VectorXd& x) const
{
if (ideq_w() != no_equation) return x(ideq_w());
else return 1.0;
}
void set_ionic_strength(const Eigen::VectorXd& x);
void compute_log_gamma(const Eigen::VectorXd& x);
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 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:
void 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