Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F66705862
reduced_system.cpp
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, 12:38
Size
11 KB
Mime Type
text/x-c
Expires
Fri, Jun 14, 12:38 (2 d)
Engine
blob
Format
Raw Data
Handle
18240783
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reduced_system.cpp
View Options
/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include <iostream>
#include "reduced_system.hpp"
#include "utils/log.hpp"
namespace
specmicp
{
inline
double
pow10
(
double
x
)
{
return
std
::
pow
(
10.0
,
x
);
}
const
double
Adebye
=
0.5092
;
const
double
Bdebye
=
0.3283
;
void
ReducedSystem
::
number_eq
(
bool
water_equation
)
{
int
neq
=
0
;
m_ideq
.
reserve
(
m_thermo
->
nb_component
()
+
m_thermo
->
nb_mineral
());
if
(
water_equation
)
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
else
{
m_ideq
.
push_back
(
no_equation
);
}
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
// Remve components that does not exist
if
(
m_tot_conc
(
i
)
==
0
and
i
!=
1
)
{
INFO
<<
"Component "
<<
m_thermo
->
label_aqueous
(
i
)
<<
"has total concentration equal to zero, we desactivate it"
<<
std
::
endl
;
m_nonactive_component
.
push_back
(
i
);
m_ideq
.
push_back
(
no_equation
);
continue
;
}
else
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
}
m_nb_free_vars
=
neq
;
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
bool
can_precipitate
=
true
;
// Remove minerals that cannot precipitate
for
(
auto
it
=
m_nonactive_component
.
begin
();
it
!=
m_nonactive_component
.
end
();
++
it
)
{
if
(
m_thermo
->
nu_min
(
m
,
*
it
)
!=
0
)
{
can_precipitate
=
false
;
break
;
// this is not a mineral that can precipitate
}
}
if
(
can_precipitate
)
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
else
{
m_ideq
.
push_back
(
no_equation
);
}
}
m_nb_tot_vars
=
neq
;
m_nb_compl_vars
=
m_nb_tot_vars
-
m_nb_free_vars
;
// aqueous species
m_active_aqueous
.
reserve
(
m_thermo
->
nb_aqueous
());
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
bool
can_exist
=
true
;
for
(
auto
it
=
m_nonactive_component
.
begin
();
it
!=
m_nonactive_component
.
end
();
++
it
)
{
if
(
m_thermo
->
nu_aq
(
j
,
*
it
)
!=
0
)
{
can_exist
=
false
;
}
}
m_active_aqueous
.
push_back
(
can_exist
);
}
}
double
ReducedSystem
::
residual_water
(
const
Eigen
::
VectorXd
&
x
)
const
{
double
res
=
m_tot_conc
(
0
)
-
mass_water
(
x
)
/
molar_mass_water
;
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
res
-=
mass_water
(
x
)
*
m_thermo
->
nu_aq
(
j
,
0
)
*
m_secondary_conc
(
j
);
}
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
res
-=
m_thermo
->
nu_min
(
m
,
0
)
*
x
(
ideq_min
(
m
));
}
return
res
;
}
double
ReducedSystem
::
residual_component
(
const
Eigen
::
VectorXd
&
x
,
int
i
)
const
{
const
int
idp
=
ideq_paq
(
i
);
double
res
=
m_tot_conc
(
i
)
-
mass_water
(
x
)
*
pow10
(
x
(
idp
));
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
not
m_active_aqueous
[
j
])
continue
;
res
-=
mass_water
(
x
)
*
m_thermo
->
nu_aq
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
res
-=
m_thermo
->
nu_min
(
m
,
i
)
*
x
(
ideq_min
(
m
));
}
return
res
;
}
double
ReducedSystem
::
residual_mineral
(
const
Eigen
::
VectorXd
&
x
,
int
m
)
const
{
double
res
=
m_thermo
->
logkr_min
(
m
);
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
m_thermo
->
nu_min
(
m
,
i
)
!=
0
)
res
-=
m_thermo
->
nu_min
(
m
,
i
)
*
(
x
(
ideq_paq
(
i
))
+
m_loggamma
(
i
));
}
return
res
;
}
void
ReducedSystem
::
get_residuals
(
const
Eigen
::
VectorXd
&
x
,
Eigen
::
VectorXd
&
residual
)
{
set_secondary_concentration
(
x
);
if
(
ideq_w
()
!=
no_equation
)
residual
(
ideq_w
())
=
residual_water
(
x
);
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
!=
no_equation
)
residual
(
ideq_paq
(
i
))
=
residual_component
(
x
,
i
);
}
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
if
(
ideq_min
(
m
)
!=
no_equation
)
residual
(
ideq_min
(
m
))
=
residual_mineral
(
x
,
m
);
}
}
//non-optimized Finite difference, for test only
//void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
// Eigen::MatrixXd& jacobian)
//{
// const int neq = total_variables();
// Eigen::VectorXd res(total_variables());
// Eigen::VectorXd perturbed_res(total_variables());
// get_residuals(x, res);
// for (int j=0; j<neq; ++j)
// {
// double h = 1e-8*std::abs(x(j));
// //h = std::copysign(h, x(j));
// if (h==0) h = 1e-8;
// double tmp = x(j);
// x(j) += h;
// h = x(j) - tmp;
// get_residuals(x, perturbed_res);
// for (int i=0; i<neq; ++i)
// {
// jacobian(i, j) = (perturbed_res(i) - res(i))/h;
// }
// x(j) = tmp;
// }
// std::cout << jacobian << std::endl;
// return;
//}
void
ReducedSystem
::
get_jacobian
(
Eigen
::
VectorXd
&
x
,
Eigen
::
MatrixXd
&
jacobian
)
{
const
double
log10
=
std
::
log
(
10.0
);
jacobian
.
setZero
();
// water
if
(
ideq_w
()
!=
no_equation
)
{
double
tmp
=
-
1
/
molar_mass_water
;
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
m_thermo
->
nu_aq
(
j
,
0
)
==
0
)
continue
;
tmp
-=
m_thermo
->
nu_aq
(
j
,
0
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
0
)
=
tmp
;
for
(
int
k
=
1
;
k
<
m_thermo
->
nb_component
();
++
k
)
{
if
(
ideq_paq
(
k
)
==
no_equation
)
continue
;
double
tmp
=
0
;
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
tmp
-=
m_thermo
->
nu_aq
(
j
,
0
)
*
m_thermo
->
nu_aq
(
j
,
k
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
ideq_paq
(
k
))
=
x
(
ideq_w
())
*
log10
*
tmp
;
}
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
jacobian
(
0
,
ideq_min
(
m
))
=
-
m_thermo
->
nu_min
(
m
,
0
);
}
}
// aqueous component
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
continue
;
const
int
idp
=
ideq_paq
(
i
);
// aqueous components
for
(
int
k
=
1
;
k
<
m_thermo
->
nb_component
();
++
k
)
{
if
(
ideq_paq
(
k
)
==
no_equation
)
continue
;
double
tmp_iip
=
0
;
if
(
k
==
i
)
tmp_iip
-=
pow10
(
x
(
ideq_paq
(
i
)))
*
log10
;
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
tmp_iip
-=
log10
*
m_thermo
->
nu_aq
(
j
,
i
)
*
m_thermo
->
nu_aq
(
j
,
k
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_paq
(
k
))
=
mass_water
(
x
)
*
tmp_iip
;
}
// minerals
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
jacobian
(
idp
,
ideq_min
(
m
))
=
-
m_thermo
->
nu_min
(
m
,
i
);
}
if
(
ideq_w
()
!=
no_equation
)
// Water
{
double
tmp_iw
=
-
std
::
pow
(
10.0
,
x
(
idp
));
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
m_thermo
->
nu_aq
(
j
,
i
)
==
0
)
continue
;
tmp_iw
-=
m_thermo
->
nu_aq
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_w
())
=
tmp_iw
;
}
}
// mineral equilibrium
for
(
int
m
=
0
;
m
<
m_thermo
->
nb_mineral
();
++
m
)
{
const
int
idm
=
ideq_min
(
m
);
if
(
idm
==
no_equation
)
continue
;
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
continue
;
jacobian
(
idm
,
ideq_paq
(
i
))
=
-
m_thermo
->
nu_min
(
m
,
i
);
}
}
//std::cout << jacobian << std::endl;
}
void
ReducedSystem
::
set_secondary_concentration
(
const
Eigen
::
VectorXd
&
x
)
{
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
m_active_aqueous
[
j
]
==
false
)
{
m_secondary_conc
(
j
)
=
0
;
continue
;
}
double
conc
=
-
m_thermo
->
logkr_aq
(
j
)
-
m_loggamma
(
j
+
m_thermo
->
nb_component
());
for
(
int
k
=
1
;
k
<
m_thermo
->
nb_component
();
++
k
)
{
if
(
m_thermo
->
nu_aq
(
j
,
k
)
==
0
)
continue
;
conc
+=
m_thermo
->
nu_aq
(
j
,
k
)
*
(
x
(
ideq_paq
(
k
))
+
m_loggamma
(
k
));
}
m_secondary_conc
(
j
)
=
pow10
(
conc
);
}
}
void
ReducedSystem
::
set_ionic_strength
(
const
Eigen
::
VectorXd
&
x
)
{
double
ionic
=
0
;
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
or
m_thermo
->
charge_paq
(
i
)
==
0
)
continue
;
ionic
+=
pow10
(
x
(
ideq_paq
(
i
)))
*
std
::
pow
(
m_thermo
->
charge_paq
(
i
),
2
);
}
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
not
m_active_aqueous
[
j
]
or
m_thermo
->
charge_saq
(
j
)
==
0
)
continue
;
ionic
+=
m_secondary_conc
(
j
)
*
std
::
pow
(
m_thermo
->
charge_saq
(
j
),
2
);
}
m_ionic_strength
=
ionic
;
}
void
ReducedSystem
::
compute_log_gamma
(
const
Eigen
::
VectorXd
&
x
)
{
set_ionic_strength
(
x
);
const
double
sqrti
=
std
::
sqrt
(
m_ionic_strength
);
for
(
int
i
=
0
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
{
m_loggamma
(
i
)
=
0
;
continue
;
}
double
tmp
=
0
;
if
(
m_thermo
->
charge
(
i
)
!=
0
)
{
tmp
=
-
Adebye
*
std
::
pow
(
2
,
m_thermo
->
charge
(
i
))
*
sqrti
;
tmp
/=
(
1
+
m_thermo
->
ao_debye
(
i
)
*
Bdebye
*
sqrti
);
}
tmp
+=
m_thermo
->
bdot_debye
(
i
)
*
m_ionic_strength
;
m_loggamma
(
i
)
=
tmp
;
}
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
not
m_active_aqueous
[
j
])
{
m_loggamma
(
m_thermo
->
nb_component
()
+
j
)
=
0
;
continue
;
}
double
tmp
=
0
;
if
(
m_thermo
->
charge_saq
(
j
)
!=
0
)
{
tmp
=
-
Adebye
*
std
::
pow
(
2
,
m_thermo
->
charge_saq
(
j
))
*
sqrti
;
tmp
/=
(
1
+
m_thermo
->
ao_debye_saq
(
j
)
*
Bdebye
*
sqrti
);
}
tmp
+=
m_thermo
->
bdot_debye_saq
(
j
)
*
m_ionic_strength
;
m_loggamma
(
m_thermo
->
nb_component
()
+
j
)
=
tmp
;
}
}
void
ReducedSystem
::
aqueous_tot_conc
(
const
Eigen
::
VectorXd
&
x
,
Eigen
::
VectorXd
&
aq_totconc
)
{
assert
(
aq_totconc
.
rows
()
==
m_thermo
->
nb_component
());
set_secondary_concentration
(
x
);
// just to be sure
for
(
int
i
=
1
;
i
<
m_thermo
->
nb_component
();
++
i
)
{
if
(
ideq_paq
(
i
)
!=
no_equation
)
{
double
tmp
=
pow10
(
x
(
ideq_paq
(
i
)));
for
(
int
j
=
0
;
j
<
m_thermo
->
nb_aqueous
();
++
j
)
{
if
(
m_thermo
->
nu_aq
(
j
,
i
)
==
0
)
continue
;
tmp
+=
m_thermo
->
nu_aq
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
aq_totconc
(
i
)
=
tmp
;
}
else
{
aq_totconc
(
i
)
=
0
;
}
}
}
double
ReducedSystem
::
max_lambda
(
const
Eigen
::
VectorXd
&
x
,
const
Eigen
::
VectorXd
&
update
)
{
//return 1.0/std::max(1.0,(-update.block(0,0,5,1).cwiseQuotient(0.95*x.block(0,0,5,1))).maxCoeff());
if
(
ideq_w
()
!=
no_equation
)
{
// std::cout << "max lambda : " << 1.0/std::max(1.0, -update(0)/(0.95*x(0))) << std::endl;
return
1.0
/
std
::
max
(
1.0
,
-
update
(
0
)
/
(
0.5
*
x
(
0
)));
}
else
{
return
1.0
;
}
}
}
// end namespace specmicp
Event Timeline
Log In to Comment