Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F66698632
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
Wed, Jun 12, 11:29
Size
9 KB
Mime Type
text/x-c++
Expires
Fri, Jun 14, 11:29 (2 d)
Engine
blob
Format
Raw Data
Handle
18157142
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Event Timeline
Log In to Comment