Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/src/specmicp/extended_system.hpp
7faccfcacb59master

Last Change
Raw File

extended_system.hpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_EXTENDEDSYSTEM_HPP
#define SPECMICP_EXTENDEDSYSTEM_HPP
//! \file extended_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 {
class SpecProg : public micpsolver::MiCPProg<SpecProg>
{
public:
const int no_equation = -1;
SpecProg(std::shared_ptr<ThermoData> ptrthermo,
Eigen::VectorXd totconc,
bool water_eq=true):
m_thermo(ptrthermo),
m_tot_conc(totconc),
m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous())
{
number_eq(water_eq);
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 {return m_ideq[i];}
int ideq_saq(int j) const {return m_ideq[j+m_thermo->nb_component()];}
int ideq_min(int m) const {return m_ideq[m+m_thermo->nb_component()+m_thermo->nb_aqueous()];}
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& _, double __) {
//compute_log_gamma(x);
not_in_linesearch = true;
}
void aqueous_tot_conc(const Eigen::VectorXd& x,
Eigen::VectorXd& aq_totconc);
private:
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;
double m_ionic_strength;
Eigen::VectorXd m_loggamma;
bool not_in_linesearch;
};
} // end namespace specmicp
#endif // SPECMICP_EXTENDEDSYSTEM_HPP
c4science · Help