Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65730166
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 5, 20:18
Size
15 KB
Mime Type
text/x-c
Expires
Fri, Jun 7, 20:18 (2 d)
Engine
blob
Format
Raw Data
Handle
18118674
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reduced_system.cpp
View Options
/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.cpp
- 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.
---------------------------------------------------------*/
#include <iostream>
#include "reduced_system.hpp"
#include "utils/log.hpp"
#include "equilibrium_data.hpp"
#include "physics/constants.hpp"
#include "physics/laws.hpp"
namespace
specmicp
{
inline
double
pow10
(
double
x
)
{
return
std
::
pow
(
10.0
,
x
);
}
ReducedSystem
::
ReducedSystem
(
std
::
shared_ptr
<
database
::
DataContainer
>
ptrdata
,
const
Eigen
::
VectorXd
&
totconc
)
:
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
ptrdata
->
nb_aqueous
),
m_loggamma
(
ptrdata
->
nb_component
+
ptrdata
->
nb_aqueous
)
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
(
true
);
m_secondary_conc
.
setZero
();
m_loggamma
.
setZero
();
}
ReducedSystem
::
ReducedSystem
(
std
::
shared_ptr
<
database
::
DataContainer
>
ptrdata
,
const
Eigen
::
VectorXd
&
totconc
,
const
SimulationBox
&
simul_box
)
:
m_simulbox
(
simul_box
),
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
ptrdata
->
nb_aqueous
),
m_loggamma
(
ptrdata
->
nb_component
+
ptrdata
->
nb_aqueous
)
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
(
true
);
m_secondary_conc
.
setZero
();
m_loggamma
.
setZero
();
}
ReducedSystem
::
ReducedSystem
(
std
::
shared_ptr
<
database
::
DataContainer
>
ptrdata
,
const
Eigen
::
VectorXd
&
totconc
,
const
SimulationBox
&
simul_box
,
const
EquilibriumState
&
previous_solution
)
:
m_simulbox
(
simul_box
),
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
previous_solution
.
molality_secondary
()),
m_loggamma
(
previous_solution
.
activity_coefficient
())
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
(
true
);
}
void
ReducedSystem
::
number_eq
(
bool
water_equation
)
{
int
neq
=
0
;
m_ideq
.
reserve
(
m_data
->
nb_component
+
m_data
->
nb_mineral
);
if
(
water_equation
)
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
else
{
m_ideq
.
push_back
(
no_equation
);
}
for
(
int
i
=
1
;
i
<
m_data
->
nb_component
;
++
i
)
{
// Remove components that does not exist
if
(
m_tot_conc
(
i
)
==
0
and
i
!=
1
)
{
INFO
<<
"Component "
<<
m_data
->
labels_basis
[
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_data
->
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_data
->
nu_mineral
(
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_data
->
nb_aqueous
);
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
bool
can_exist
=
true
;
for
(
auto
it
=
m_nonactive_component
.
begin
();
it
!=
m_nonactive_component
.
end
();
++
it
)
{
if
(
m_data
->
nu_aqueous
(
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
)
/
m_data
->
molar_mass_basis_si
(
0
);
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
m_data
->
nu_aqueous
(
j
,
0
)
==
0
)
continue
;
res
-=
mass_water
(
x
)
*
m_data
->
nu_aqueous
(
j
,
0
)
*
m_secondary_conc
(
j
);
}
for
(
int
m
=
0
;
m
<
m_data
->
nb_mineral
;
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
or
m_data
->
nu_mineral
(
m
,
0
)
==
0
)
continue
;
res
-=
m_data
->
nu_mineral
(
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_data
->
nb_aqueous
;
++
j
)
{
if
(
not
m_active_aqueous
[
j
])
continue
;
res
-=
mass_water
(
x
)
*
m_data
->
nu_aqueous
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
for
(
int
m
=
0
;
m
<
m_data
->
nb_mineral
;
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
res
-=
m_data
->
nu_mineral
(
m
,
i
)
*
x
(
ideq_min
(
m
));
}
return
res
;
}
double
ReducedSystem
::
residual_mineral
(
const
Eigen
::
VectorXd
&
x
,
int
m
)
const
{
double
res
=
m_data
->
logk_mineral
(
m
);
for
(
int
i
=
1
;
i
<
m_data
->
nb_component
;
++
i
)
{
if
(
m_data
->
nu_mineral
(
m
,
i
)
!=
0
)
res
-=
m_data
->
nu_mineral
(
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_data
->
nb_component
;
++
i
)
{
if
(
ideq_paq
(
i
)
!=
no_equation
)
residual
(
ideq_paq
(
i
))
=
residual_component
(
x
,
i
);
}
for
(
int
m
=
0
;
m
<
m_data
->
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
.
resize
(
total_variables
(),
total_variables
());
jacobian
.
setZero
();
// water
if
(
ideq_w
()
!=
no_equation
)
{
double
tmp
=
-
1
/
m_data
->
molar_mass_basis_si
(
0
);
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
m_data
->
nu_aqueous
(
j
,
0
)
==
0
)
continue
;
tmp
-=
m_data
->
nu_aqueous
(
j
,
0
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
0
)
=
tmp
;
for
(
int
k
=
1
;
k
<
m_data
->
nb_component
;
++
k
)
{
if
(
ideq_paq
(
k
)
==
no_equation
)
continue
;
double
tmp
=
0
;
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
tmp
-=
m_data
->
nu_aqueous
(
j
,
0
)
*
m_data
->
nu_aqueous
(
j
,
k
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
ideq_paq
(
k
))
=
x
(
ideq_w
())
*
log10
*
tmp
;
}
for
(
int
m
=
0
;
m
<
m_data
->
nb_mineral
;
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
jacobian
(
0
,
ideq_min
(
m
))
=
-
m_data
->
nu_mineral
(
m
,
0
);
}
}
// aqueous component
for
(
int
i
=
1
;
i
<
m_data
->
nb_component
;
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
continue
;
const
int
idp
=
ideq_paq
(
i
);
// aqueous components
for
(
int
k
=
1
;
k
<
m_data
->
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_data
->
nb_aqueous
;
++
j
)
{
tmp_iip
-=
log10
*
m_data
->
nu_aqueous
(
j
,
i
)
*
m_data
->
nu_aqueous
(
j
,
k
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_paq
(
k
))
=
mass_water
(
x
)
*
tmp_iip
;
}
// minerals
for
(
int
m
=
0
;
m
<
m_data
->
nb_mineral
;
++
m
)
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
jacobian
(
idp
,
ideq_min
(
m
))
=
-
m_data
->
nu_mineral
(
m
,
i
);
}
if
(
ideq_w
()
!=
no_equation
)
// Water
{
double
tmp_iw
=
-
std
::
pow
(
10.0
,
x
(
idp
));
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
m_data
->
nu_aqueous
(
j
,
i
)
==
0
)
continue
;
tmp_iw
-=
m_data
->
nu_aqueous
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_w
())
=
tmp_iw
;
}
}
// mineral equilibrium
for
(
int
m
=
0
;
m
<
m_data
->
nb_mineral
;
++
m
)
{
const
int
idm
=
ideq_min
(
m
);
if
(
idm
==
no_equation
)
continue
;
for
(
int
i
=
1
;
i
<
m_data
->
nb_component
;
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
continue
;
jacobian
(
idm
,
ideq_paq
(
i
))
=
-
m_data
->
nu_mineral
(
m
,
i
);
}
}
}
void
ReducedSystem
::
set_secondary_concentration
(
const
Eigen
::
VectorXd
&
x
)
{
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
m_active_aqueous
[
j
]
==
false
)
{
m_secondary_conc
(
j
)
=
0
;
continue
;
}
double
conc
=
-
m_data
->
logk_aqueous
(
j
)
-
m_loggamma
(
j
+
m_data
->
nb_component
);
for
(
int
k
=
1
;
k
<
m_data
->
nb_component
;
++
k
)
{
if
(
m_data
->
nu_aqueous
(
j
,
k
)
==
0
)
continue
;
conc
+=
m_data
->
nu_aqueous
(
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_data
->
nb_component
;
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
or
m_data
->
charge_component
(
i
)
==
0
)
continue
;
ionic
+=
pow10
(
x
(
ideq_paq
(
i
)))
*
std
::
pow
(
m_data
->
charge_component
(
i
),
2
);
}
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
not
m_active_aqueous
[
j
]
or
m_data
->
charge_aqueous
(
j
)
==
0
)
continue
;
ionic
+=
m_secondary_conc
(
j
)
*
std
::
pow
(
m_data
->
charge_aqueous
(
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_data
->
nb_component
;
++
i
)
{
if
(
ideq_paq
(
i
)
==
no_equation
)
{
m_loggamma
(
i
)
=
0
;
continue
;
}
m_loggamma
(
i
)
=
laws
::
extended_debye_huckel
(
m_ionic_strength
,
sqrti
,
m_data
->
charge_component
(
i
),
m_data
->
a_debye_component
(
i
),
m_data
->
b_debye_component
(
i
));
}
for
(
int
j
=
0
;
j
<
m_data
->
nb_aqueous
;
++
j
)
{
if
(
not
m_active_aqueous
[
j
])
{
m_loggamma
(
m_data
->
nb_component
+
j
)
=
0
;
continue
;
}
m_loggamma
(
m_data
->
nb_component
+
j
)
=
laws
::
extended_debye_huckel
(
m_ionic_strength
,
sqrti
,
m_data
->
charge_aqueous
(
j
),
m_data
->
a_debye_aqueous
(
j
),
m_data
->
b_debye_aqueous
(
j
));
}
}
bool
ReducedSystem
::
hook_start_iteration
(
const
Eigen
::
VectorXd
&
x
,
double
norm_residual
)
{
not_in_linesearch
=
true
;
double
previous_norm
=
m_loggamma
.
norm
();
bool
may_have_converged
=
false
;
if
(
norm_residual
<
nb_free_variables
())
{
// Use fix point iterations for non-ideality
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
set_secondary_concentration
(
x
);
compute_log_gamma
(
x
);
if
(
std
::
abs
(
previous_norm
-
m_loggamma
.
norm
())
<
1e-6
)
{
may_have_converged
=
true
;
break
;}
previous_norm
=
m_loggamma
.
norm
();
}
}
return
may_have_converged
;
}
double
ReducedSystem
::
max_lambda
(
const
Eigen
::
VectorXd
&
x
,
const
Eigen
::
VectorXd
&
update
)
{
if
(
ideq_w
()
!=
no_equation
)
{
return
1.0
/
std
::
max
(
1.0
,
-
update
(
0
)
/
(
0.9
*
x
(
0
)));
}
else
{
return
1.0
;
}
}
void
ReducedSystem
::
reasonable_starting_guess
(
Eigen
::
VectorXd
&
x
,
bool
for_min
)
{
x
.
resize
(
m_data
->
nb_component
+
m_data
->
nb_mineral
);
x
(
0
)
=
2
*
m_tot_conc
(
0
)
*
m_data
->
molar_mass_basis_si
(
0
);
for
(
int
i
=
1
;
i
<
m_data
->
nb_component
;
++
i
)
{
x
(
i
)
=
-
4.0
;
}
if
(
for_min
)
x
.
block
(
m_data
->
nb_component
,
0
,
m_data
->
nb_mineral
,
1
).
setZero
();
else
x
.
block
(
m_data
->
nb_component
,
0
,
m_data
->
nb_mineral
,
1
).
setZero
();
m_loggamma
.
setZero
();
m_secondary_conc
.
setZero
();
}
void
ReducedSystem
::
reasonable_restarting_guess
(
Eigen
::
VectorXd
&
x
)
{
x
(
0
)
=
m_tot_conc
(
0
)
*
m_data
->
molar_mass_basis_si
(
0
);
for
(
int
i
=
m_data
->
nb_component
;
i
<
m_data
->
nb_component
+
m_data
->
nb_mineral
;
++
i
)
{
x
(
i
)
=
0.0
;
}
m_loggamma
.
setZero
();
m_secondary_conc
.
setZero
();
}
EquilibriumState
ReducedSystem
::
get_solution
(
const
Eigen
::
VectorXd
&
xtot
,
const
Eigen
::
VectorXd
&
x
)
{
double
previous_norm
=
m_loggamma
.
norm
();
set_secondary_concentration
(
x
);
compute_log_gamma
(
x
);
if
(
std
::
abs
(
previous_norm
-
m_loggamma
.
norm
())
>
1e-6
)
{
WARNING
<<
"Activity coefficient have not converged !"
<<
std
::
endl
<<
"output can not be trusted
\n
Difference "
+
std
::
to_string
(
std
::
abs
(
previous_norm
-
m_loggamma
.
norm
()));
}
return
EquilibriumState
(
xtot
,
m_secondary_conc
,
m_loggamma
,
m_ionic_strength
,
m_data
);
}
}
// end namespace specmicp
Event Timeline
Log In to Comment