Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62184095
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
Sat, May 11, 11:34
Size
3 KB
Mime Type
text/x-c++
Expires
Mon, May 13, 11:34 (2 d)
Engine
blob
Format
Raw Data
Handle
17615964
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