Homec4science
Diffusion SpecMiCP / ReactMiCP (oxydoreduc) Browse

specmicp/src/specmicp/reduced_system.hpp
d442c9ff39a3oxydoreduc

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 "common.hpp"
#include "database.hpp"
#include "micpsolver/micpprog.hpp"
#include "simulation_box.hpp"
#include "reduced_system_structs.hpp"
#include "utils/options_handler.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 OptionsHandler<ReducedSystemOptions>
{
public:
//! \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(
RawDatabasePtr ptrdata,
const Vector& totconc,
const ReducedSystemOptions& options
);
//! \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(
RawDatabasePtr ptrdata,
const Vector& totconc,
const SimulationBox& simul_box,
const ReducedSystemOptions& options
);
ReducedSystem(
RawDatabasePtr ptrdata,
const Vector& totconc,
const SimulationBox& simul_box,
const EquilibriumState& previous_solution,
const ReducedSystemOptions& options
);
// Variables
// ==========
//! \brief Return the total number of variables
index_t total_variables() const {return m_nb_tot_vars;}
//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
index_t nb_free_variables() const {return m_nb_free_vars;}
//! \brief Return the number of variables subjected to the complementarity conditions
index_t nb_complementarity_variables() const {return m_nb_compl_vars;}
// The linear system
// ==================
//! \brief Return the residual for the water conservation equation
scalar_t residual_water(const Vector& 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 Vector& x, index_t i) const;
//! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$
double residual_mineral(const Vector& x, index_t m) const;
//! \brief Charge conservation
double residual_charge_conservation(const Vector& x) const;
//! \brief Return the residuals
void get_residuals(const Vector& x, Vector& residual);
//! \brief Return the jacobian
void get_jacobian(Vector& x,
Matrix& jacobian);
//! \brief Return the equation id
index_t ideq(index_t i) const {return m_ideq[i];}
//! \brief Return the equation number for conservation of water
index_t ideq_w() const {return m_ideq[0];}
//! \brief Return the equation number of component 'i'
index_t ideq_paq(index_t i) const {
specmicp_assert(i < m_data->nb_component);
return m_ideq[i];
}
//! \brief Return the equation number of mineral 'm'
index_t ideq_min(index_t m) const {
specmicp_assert(m < m_data->nb_mineral);
return m_ideq[m+m_data->nb_component];
}
//! \brief Return the mass of water
scalar_t mass_water(const Vector& x) const
{
if (ideq_w() != no_equation) return x(ideq_w());
else return 1.0; //FIXME
}
scalar_t component_concentration(const Vector& x, index_t component) const
{
specmicp_assert(ideq_paq(component) != no_equation);
return std::pow(10.0, x(ideq_paq(component)));
}
//! \brief Compute the ionic strength
void set_ionic_strength(const Vector& x);
//! \brief Compute the activity coefficients
void compute_log_gamma(const Vector& x);
//! \brief Compute the activity model, called by the solver at the beginning of an iteration
bool hook_start_iteration(const Vector &x, scalar_t norm_residual);
//! \brief Return a scale factor to avoid negative mass during Newton's iteration
double max_lambda(const Vector &x, const Vector &update);
//! \brief Return the component that does not exist in the solution (inactive dof)
const std::vector<index_t>& get_non_active_component() {return m_nonactive_component;}
//! \brief A reasonable (well... maybe...) starting guess
void reasonable_starting_guess(Vector& x, bool for_min=false);
//! \brief A reasonable (maybe...) restarting guess
void reasonable_restarting_guess(Vector& x);
//! Return the equilibrium state of the system, the 'Solution' of the speciation problem
EquilibriumState get_solution(const Vector& xtot, const Vector& x);
//! \brief Return the total concentration of 'component'
scalar_t total_concentration(index_t component) const {return m_tot_conc(component);}
//! \brief Reduce the mineral system if one cannot dissolve
void reduce_mineral_problem(Vector& true_var);
//! \brief Reset the mineral system to the true problem
void reset_mineral_system(Vector& true_var, Vector& output_var);
//! \brief set the charge keeper to 'component'
//!
//! Set component to 'no_species' to solve charge conservation implicitely.
//! In this case, the total concentration must conserve the charge.
void set_charge_keeper(index_t component)
{
specmicp_assert(component >= 1 and component < m_data->nb_component);
m_charge_keeper = component;
}
//! Return true if 'component' is thecharge keeper
bool is_charge_keeper(index_t component)
{
return (component == m_charge_keeper);
}
private:
void jacobian_water(Vector& x, Matrix& jacobian);
void jacobian_aqueous_components(Vector& x, Matrix& jacobian);
void jacobian_minerals(Vector& x, Matrix& jacobian);
// For the normal algorithm
void set_secondary_concentration(const Vector& x);
// For the tot aq_conc pb
//void unsafe_set_secondary_concentration(const Eigen::VectorXd& x);
void number_eq();
SimulationBox m_simulbox;
RawDatabasePtr m_data;
Vector m_tot_conc;
index_t m_nb_tot_vars;
index_t m_nb_free_vars;
index_t m_nb_compl_vars;
std::vector<index_t> m_ideq;
Vector m_secondary_conc;
Vector m_gas_pressure;
scalar_t m_ionic_strength;
Vector m_loggamma;
index_t m_charge_keeper;
bool not_in_linesearch;
std::vector<index_t> m_nonactive_component;
std::vector<bool> m_active_aqueous;
std::vector<bool> m_active_gas;
};
} // end namespace specmicp
#endif // SPECMICP_REDUCEDSYSTEM_HPP
c4science · Help