Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F67209861
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
Thu, Jun 20, 20:19
Size
7 KB
Mime Type
text/x-c++
Expires
Sat, Jun 22, 20:19 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18350961
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 "specmicp/reduced_system_solver.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
;
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
<
int
>
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
/
specmicp
::
molar_mass_water
,
1e-7
,
1e-6
,
2e-6
;
m_totconc0
<<
1.0
/
specmicp
::
molar_mass_water
,
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
);
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
;
};
int
main
()
{
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
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;
return
EXIT_SUCCESS
;
}
Event Timeline
Log In to Comment