Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/reduced_system.hpp
be84b355a07dreactivetransport

Last Change
Raw File

reduced_system.hpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_REDUCEDSYSTEM_HPP
#define SPECMICP_REDUCEDSYSTEM_HPP
//! \file reduced_system.hpp The MiCP program to solve speciation
#include <memory>
#include "thermodata.hpp"
#include "kineticsdata.hpp"
#include "micpsolver/micpprog.hpp"
#include "simulation_box.hpp"
//! \namespace specmicp main namespace for the SpecMiCP solver
namespace specmicp {
class EquilibriumState;
//! \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
ReducedSystem(std::shared_ptr<ThermoData> ptrthermo,
const Eigen::VectorXd& totconc);
//! \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
ReducedSystem(std::shared_ptr<ThermoData> ptrthermo,
const Eigen::VectorXd& totconc,
const SimulationBox& simul_box);
ReducedSystem(std::shared_ptr<ThermoData> ptrthermo,
const Eigen::VectorXd& totconc,
const SimulationBox& simul_box,
const EquilibriumState& previous_solution);
// 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
bool hook_start_iteration(const Eigen::VectorXd& x, double norm_residual);
//! \brief Return a scale factor to avoid negative mass during Newton's iteration
double max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update);
//! \brief Return the component that does not exist in the solution (inactive dof)
const std::vector<int>& get_non_active_component() {return m_nonactive_component;}
//! \brief A reasonable (well... maybe...) starting guess
void reasonable_starting_guess(Eigen::VectorXd& x);
//! \brief A reasonable (maybe...) restarting guess
void reasonable_restarting_guess(Eigen::VectorXd& x);
//! Return the equilibrium state of the system, the 'Solution' of the speciation problem
EquilibriumState get_solution(const Eigen::VectorXd& xtot, const 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);
SimulationBox m_simulbox;
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;
Eigen::VectorXd m_gas_pressure;
double m_ionic_strength;
Eigen::VectorXd m_loggamma;
bool not_in_linesearch;
std::vector<int> m_nonactive_component;
std::vector<bool> m_active_aqueous;
std::vector<bool> m_active_gas;
};
} // end namespace specmicp
#endif // SPECMICP_REDUCEDSYSTEM_HPP
c4science · Help