Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74572078
reduced_system.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Jul 28, 12:38
Size
3 KB
Mime Type
text/x-c++
Expires
Tue, Jul 30, 12:38 (2 d)
Engine
blob
Format
Raw Data
Handle
19410501
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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:
// 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
Event Timeline
Log In to Comment