Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71248967
kinetics.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, Jul 10, 13:54
Size
13 KB
Mime Type
text/x-c++
Expires
Fri, Jul 12, 13:54 (2 d)
Engine
blob
Format
Raw Data
Handle
18926976
Attached To
rSPECMICP SpecMiCP / ReactMiCP
kinetics.cpp
View Options
/*-------------------------------------------------------
- Module : tests/
- File : kinetics.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 "common.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include "database.hpp"
//#include "database/database.hpp"
#include <iostream>
#include <iomanip>
#include "utils/log.hpp"
#include "odeint/runge_kutta_step.hpp"
#include "specmicp/equilibrium_data.hpp"
const
double
Mc3s
=
228.3
;
const
double
Vmc3s
=
73e-6
;
using
specmicp
::
index_t
;
using
specmicp
::
scalar_t
;
class
RateComputer
{
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
public
:
static
constexpr
double
m_c3s0
=
500
;
const
double
R0
=
2.5e-6
;
RateComputer
()
:
m_database
(
"../data/cemdata_specmicp.js"
)
{
std
::
shared_ptr
<
specmicp
::
database
::
DataContainer
>
data
=
m_database
.
get_database
();
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
}
});
m_database
.
swap_components
(
swapping
);
std
::
vector
<
index_t
>
to_keep
=
{
0
,
1
,
m_database
.
safe_label_to_id
(
"Ca[2+]"
,
m_database
.
get_database
()
->
labels_basis
),
m_database
.
safe_label_to_id
(
"SiO(OH)3[-]"
,
m_database
.
get_database
()
->
labels_basis
)
};
m_database
.
keep_only_components
(
to_keep
);
m_equilvar
.
resize
(
data
->
nb_component
+
data
->
nb_mineral
);
m_equilvar
(
0
)
=
1.0
;
m_equilvar
.
block
(
1
,
0
,
data
->
nb_component
-
1
,
1
).
setConstant
(
-
7
);
m_equilvar
.
block
(
data
->
nb_component
,
0
,
data
->
nb_mineral
,
1
).
setConstant
(
0.
);
m_totconc
=
Eigen
::
VectorXd
(
data
->
nb_component
);
m_totconc0
=
Eigen
::
VectorXd
(
data
->
nb_component
);
m_totconc
<<
1.0
/
data
->
molar_mass_basis
(
0
),
1e-7
,
1e-6
,
2e-6
;
m_totconc0
<<
1.0
/
data
->
molar_mass_basis
(
0
),
1e-7
,
1e-6
,
2e-6
;
id_c3s
=
m_database
.
mineral_kinetic_label_to_id
(
"C3S"
);
id_ch
=
m_database
.
mineral_label_to_id
(
"Portlandite"
);
id_csht
=
m_database
.
mineral_label_to_id
(
"CSH,tobermorite"
);
id_cshj
=
m_database
.
mineral_label_to_id
(
"CSH,jennite"
);
id_ca
=
m_database
.
component_label_to_id
(
"Ca[2+]"
);
n_particle
=
3
*
Vmc3s
*
m_c3s0
/
(
Mc3s
*
4
*
M_PI
*
std
::
pow
(
R0
,
3.0
));
}
void
update_totconc
(
double
mc3s
)
{
std
::
shared_ptr
<
specmicp
::
database
::
DataContainer
>
data
=
m_database
.
get_database
();
double
xi
=
(
m_c3s0
-
mc3s
)
/
Mc3s
;
std
::
cout
<<
"xi : "
<<
xi
<<
std
::
endl
;
//std::cout <<xi*data->nu_mineral_kinetic.row(id_c3s).transpose() << std::endl;
m_totconc
.
noalias
()
=
m_totconc0
+
xi
*
data
->
nu_mineral_kinetic
.
row
(
id_c3s
).
transpose
();
}
double
compute_logiap
()
{
return
m_current_solution
.
logIAP_kinetic
(
id_c3s
);
}
double
compute_rate
(
double
mc3s
,
double
logiap
)
{
double
radius
=
std
::
pow
((
Vmc3s
*
mc3s
*
3
)
/
(
Mc3s
*
n_particle
*
4
*
M_PI
),
1.0
/
3.0
);
double
surf
=
(
3
*
Vmc3s
/
(
Mc3s
*
radius
));
std
::
cout
<<
n_particle
<<
"
\t
"
<<
radius
<<
"
\t
"
<<
surf
<<
std
::
endl
;
return
std
::
min
((
surf
*
mc3s
)
*
Mc3s
*
1e-6
*
(
2.1
*
logiap
+
35.7
),
-
1e-6
);
}
void
get_rate
(
double
_
,
double
mc3s
,
double
&
rate
)
{
//std::cout << m_totconc << std::endl;
update_totconc
(
mc3s
);
//std::cout << m_totconc << std::endl;
compute_equilibrium
();
rate
=
compute_rate
(
mc3s
,
compute_logiap
());
}
double
get_ca
()
{
return
m_current_solution
.
molality_component
(
id_ca
);
}
double
get_ch
()
{
return
m_current_solution
.
mass_mineral
(
id_ch
);
}
double
get_csht
()
{
return
m_current_solution
.
mass_mineral
(
id_csht
);
}
double
get_cshj
()
{
return
m_current_solution
.
mass_mineral
(
id_cshj
);
}
double
get_pH
()
{
return
m_current_solution
.
pH
();
}
void
compute_equilibrium
()
{
specmicp
::
ReducedSystemSolver
solver
(
m_database
.
get_database
(),
m_totconc
);
solver
.
get_options
().
solver_options
.
maxstep
=
50
;
solver
.
get_options
().
system_options
.
charge_keeper
=
1
;
specmicp
::
micpsolver
::
MiCPPerformance
perf
=
solver
.
solve
(
m_equilvar
);
if
(
perf
.
return_code
!=
specmicp
::
micpsolver
::
MiCPSolverReturnCode
::
ResidualMinimized
)
{
std
::
cout
<<
"Failed to solve the system ! Err code "
<<
(
int
)
perf
.
return_code
<<
std
::
endl
;
}
m_current_solution
=
solver
.
get_solution
(
m_equilvar
);
}
Eigen
::
VectorXd
get_vars
()
{
return
m_equilvar
;}
private
:
Eigen
::
VectorXd
m_equilvar
;
Eigen
::
VectorXd
m_totconc0
;
Eigen
::
VectorXd
m_totconc
;
specmicp
::
database
::
Database
m_database
;
specmicp
::
EquilibriumState
m_current_solution
;
int
id_c3s
;
int
id_ch
;
int
id_csht
;
int
id_cshj
;
int
id_ca
;
double
n_particle
;
};
void
by_hand
()
{
RateComputer
rate
;
double
y
=
RateComputer
::
m_c3s0
-
0.05
;
double
dydx
;
rate
.
get_rate
(
0
,
y
,
dydx
);
std
::
cout
<<
std
::
setprecision
(
6
);
specmicp
::
odeint
::
rhs_f
<
1
>
rhs
=
std
::
bind
(
std
::
mem_fn
(
&
RateComputer
::
get_rate
),
rate
,
std
::
placeholders
::
_1
,
std
::
placeholders
::
_2
,
std
::
placeholders
::
_3
);
specmicp
::
odeint
::
CashKarpStep
<
1
>
driver
=
specmicp
::
odeint
::
get_cash_karp_step
<
1
>
(
rhs
);
//specmicp::odeint::DormandPrinceStep<1> driver = specmicp::odeint::get_dormand_prince_step<1>(rhs);
specmicp
::
odeint
::
StepLength
stepl
(
1.0
);
double
x
=
0.0
;
int
cnt
=
0
;
while
(
y
>
1e-2
)
//while (cnt < 10)
{
driver
.
rk_step
(
y
,
dydx
,
x
,
stepl
);
cnt
+=
1
;
stepl
.
get_next
();
driver
.
get_rhs
(
x
,
y
,
dydx
);
rate
.
update_totconc
(
y
);
// shouldn't be needed
rate
.
compute_equilibrium
();
std
::
cout
<<
x
<<
"
\t
"
<<
y
<<
"
\t
"
<<
rate
.
compute_logiap
()
<<
"
\t
"
<<
rate
.
get_pH
()
<<
"
\t
"
<<
rate
.
get_ca
()
<<
"
\t
"
<<
rate
.
get_ch
()
<<
"
\t
"
<<
rate
.
get_csht
()
<<
"
\t
"
<<
rate
.
get_cshj
()
<<
std
::
endl
;
}
// Eigen::VectorXd equil = rate.get_vars() ;
// y+= 0.00001;
// rate.get_rate(0, y, dydx);
// Eigen::VectorXd equil2 = rate.get_vars();
// std::cout << equil2 - equil << std::endl;
//std::cout << "logiap : " << rate.compute_logiap() << " - rate : " << rate.compute_rate(rate.compute_logiap()) << std::endl;
}
/////// using the new solver
///
#include "specmicp/kinetics/kinetic_model.hpp"
#include "specmicp/kinetics/kinetic_variables.hpp"
#include "specmicp/kinetics/kinetic_system_solver.hpp"
class
C3SDissolution
:
public
specmicp
::
kinetics
::
KineticModel
{
public
:
static
constexpr
scalar_t
R0
=
2.5e-6
;
C3SDissolution
(
specmicp
::
RawDatabasePtr
data
,
scalar_t
n_c3s0
)
{
specmicp
::
Database
db
(
data
);
add_equation
(
db
.
mineral_kinetic_label_to_id
(
"C3S"
));
m_n_c3s0
=
n_c3s0
;
n_particle
=
3.0
/
(
4.0
*
M_PI
*
std
::
pow
(
R0
,
3.0
))
*
Vmc3s
*
n_c3s0
;
}
void
compute_rate
(
scalar_t
_
,
const
specmicp
::
Vector
&
y
,
specmicp
::
kinetics
::
KineticVariables
&
variables
,
specmicp
::
Vector
&
dydt
)
{
specmicp
::
Vector
y_temp
=
y
;
//specmicp::odeint::set_non_negative<Eigen::Dynamic>(y_temp);
//std::cout << "y_temp : " << y_temp << std::endl;
dydt
.
resizeLike
(
y
);
double
radius
=
R0
*
std
::
pow
((
y_temp
(
0
)
/
m_n_c3s0
),
1.0
/
3.0
);
if
(
radius
<
0.0
or
not
std
::
isfinite
(
radius
))
radius
=
0.0
;
double
surf
=
n_particle
*
4
*
M_PI
*
radius
*
radius
;
//std::cout << n_particle << "\t" << radius << "\t" << surf << std::endl;
dydt
(
0
)
=
1e-9
*
std
::
min
(
surf
*
(
2.1
*
variables
.
equilibrium_state
().
logIAP_kinetic
(
index_species
(
0
))
+
35.7
),
-
1e-6
);
}
private
:
double
m_n_c3s0
;
double
n_particle
;
};
void
using_kinetic_solver
()
{
specmicp
::
Database
the_database
(
"../data/cemdata_specmicp.js"
);
std
::
shared_ptr
<
specmicp
::
database
::
DataContainer
>
data
=
the_database
.
get_database
();
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
}
});
the_database
.
swap_components
(
swapping
);
std
::
vector
<
index_t
>
to_keep
=
{
0
,
1
,
the_database
.
safe_label_to_id
(
"Ca[2+]"
,
the_database
.
get_database
()
->
labels_basis
),
the_database
.
safe_label_to_id
(
"SiO(OH)3[-]"
,
the_database
.
get_database
()
->
labels_basis
)
};
the_database
.
keep_only_components
(
to_keep
);
specmicp
::
Vector
totconc
(
data
->
nb_component
);
totconc
(
the_database
.
component_label_to_id
(
"H2O"
))
=
1.0
/
data
->
molar_mass_basis_si
(
0
);
totconc
(
the_database
.
component_label_to_id
(
"HO[-]"
))
=
1e-7
;
totconc
(
the_database
.
component_label_to_id
(
"Ca[2+]"
))
=
1e-4
;
totconc
(
the_database
.
component_label_to_id
(
"SiO(OH)3[-]"
))
=
2e-4
;
specmicp
::
ReducedSystemSolver
first_solver
(
data
,
totconc
);
specmicp
::
Vector
minerals
(
1
);
minerals
<<
1e-3
*
(
500
/
Mc3s
);
std
::
cout
<<
"Start : "
<<
minerals
(
0
)
<<
std
::
endl
;
specmicp
::
Vector
equilvar
;
equilvar
.
resize
(
data
->
nb_component
+
data
->
nb_mineral
);
equilvar
(
0
)
=
1.0
;
equilvar
.
block
(
1
,
0
,
data
->
nb_component
-
1
,
1
).
setConstant
(
-
7
);
equilvar
.
block
(
data
->
nb_component
,
0
,
data
->
nb_mineral
,
1
).
setConstant
(
0.
);
first_solver
.
solve
(
equilvar
);
specmicp
::
Vector
second_conc
(
data
->
nb_aqueous
);
second_conc
.
setZero
();
specmicp
::
Vector
loggamma
(
data
->
nb_component
+
data
->
nb_aqueous
);
loggamma
.
setZero
();
specmicp
::
EquilibriumState
equilsolution
=
first_solver
.
get_solution
(
equilvar
);
std
::
cout
<<
equilvar
<<
std
::
endl
;
std
::
cout
<<
"Initialization finished "
<<
std
::
endl
;
std
::
shared_ptr
<
specmicp
::
kinetics
::
KineticModel
>
model
=
std
::
make_shared
<
C3SDissolution
>
(
data
,
minerals
(
0
));
specmicp
::
kinetics
::
KineticSystemSolver
solver
(
model
,
totconc
,
minerals
,
equilsolution
);
solver
.
get_options
().
speciation_options
.
solver_options
.
maxstep
=
50
;
solver
.
get_options
().
speciation_options
.
solver_options
.
use_scaling
=
true
;
solver
.
get_options
().
speciation_options
.
solver_options
.
fvectol
=
1e-8
;
solver
.
get_options
().
speciation_options
.
system_options
.
charge_keeper
=
the_database
.
component_label_to_id
(
"HO[-]"
);
solver
.
get_options
().
speciation_options
.
solver_options
.
non_monotone_linesearch
=
true
;
solver
.
get_options
().
speciation_options
.
solver_options
.
projection_min_variable
=
1e-8
;
//solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1.0;
//solver.get_options().speciation_options.conservation_water = false;
solver
.
get_options
().
rk_options
.
non_negativity
=
true
;
solver
.
solve
(
1000.0
,
100000.0
);
// solver.solve(1.0, 2000.0);
// solver.solve(1.0, 2000.0);
std
::
cout
<<
"Solution : "
<<
solver
.
variables
().
moles_mineral
(
0
)
<<
std
::
endl
;
std
::
cout
<<
"flux : "
<<
std
::
endl
<<
solver
.
variables
().
rate_components
()
<<
std
::
endl
;
solver
.
variables
().
reset_rate
();
solver
.
solve
(
solver
.
current_dt
(),
100000.0
);
// solver.solve(1.0, 2000.0);
// solver.solve(1.0, 2000.0);
std
::
cout
<<
"Solution : "
<<
solver
.
variables
().
moles_mineral
(
0
)
<<
std
::
endl
;
std
::
cout
<<
"flux : "
<<
std
::
endl
<<
solver
.
variables
().
rate_components
()
<<
std
::
endl
;
solver
.
variables
().
reset_rate
();
solver
.
solve
(
solver
.
current_dt
(),
100000.0
);
// solver.solve(1.0, 2000.0);
// solver.solve(1.0, 2000.0);
std
::
cout
<<
"Solution : "
<<
solver
.
variables
().
moles_mineral
(
0
)
<<
std
::
endl
;
std
::
cout
<<
"flux : "
<<
std
::
endl
<<
solver
.
variables
().
rate_components
()
<<
std
::
endl
;
solver
.
variables
().
reset_rate
();
solver
.
solve
(
solver
.
current_dt
(),
100000.0
);
// solver.solve(1.0, 2000.0);
// solver.solve(1.0, 2000.0);
std
::
cout
<<
"Solution : "
<<
solver
.
variables
().
moles_mineral
(
0
)
<<
std
::
endl
;
std
::
cout
<<
"flux : "
<<
std
::
endl
<<
solver
.
variables
().
rate_components
()
<<
std
::
endl
;
}
int
main
()
{
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
using_kinetic_solver
();
//by_hand();
return
EXIT_SUCCESS
;
}
Event Timeline
Log In to Comment