Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81202396
equilibrium_curve.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, Sep 4, 17:48
Size
12 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 17:48 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20529580
Attached To
rSPECMICP SpecMiCP / ReactMiCP
equilibrium_curve.cpp
View Options
#include <iostream>
#include "utils/log.hpp"
#include "reactmicp/equilibrium_curve/chemistry.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_extractor.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_coupler.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_solid_transport.hpp"
specmicp
::
Matrix
test_chemistry
()
{
specmicp
::
database
::
Database
thedatabase
(
"../data/cemdata_specmicp.js"
);
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
},
});
thedatabase
.
swap_components
(
swapping
);
thedatabase
.
remove_gas_phases
();
specmicp
::
RawDatabasePtr
raw_data
=
thedatabase
.
get_database
();
specmicp
::
Formulation
formulation
;
specmicp
::
scalar_t
mult
=
6.5e3
;
specmicp
::
scalar_t
m_c3s
=
mult
*
0.7
;
specmicp
::
scalar_t
m_c2s
=
mult
*
0.3
;
specmicp
::
scalar_t
wc
=
0.5
;
specmicp
::
scalar_t
m_water
=
wc
*
1e-3
*
(
m_c3s
*
(
3
*
56.08
+
60.08
)
+
m_c2s
*
(
2
*
56.06
+
60.08
)
);
formulation
.
mass_solution
=
m_water
;
formulation
.
amount_minerals
=
{
{
"C3S"
,
m_c3s
},
{
"C2S"
,
m_c2s
},
};
specmicp
::
Vector
total_concentrations
=
specmicp
::
Dissolver
(
raw_data
).
dissolve
(
formulation
);
specmicp
::
index_t
id_h2o
=
thedatabase
.
component_label_to_id
(
"H2O"
);
specmicp
::
index_t
id_ho
=
thedatabase
.
component_label_to_id
(
"HO[-]"
);
specmicp
::
index_t
id_ca
=
thedatabase
.
component_label_to_id
(
"Ca[2+]"
);
specmicp
::
AdimensionalSystemBC
conditions
(
total_concentrations
);
conditions
.
charge_keeper
=
id_ho
;
specmicp
::
AdimensionalSystemSolverOptions
options
;
options
.
solver_options
.
maxstep
=
10.0
;
options
.
solver_options
.
max_iter
=
100
;
options
.
solver_options
.
maxiter_maxstep
=
100
;
options
.
solver_options
.
use_crashing
=
false
;
options
.
solver_options
.
use_scaling
=
false
;
options
.
solver_options
.
factor_descent_condition
=
-
1
;
options
.
solver_options
.
factor_gradient_search_direction
=
100
;
options
.
solver_options
.
projection_min_variable
=
1e-9
;
options
.
solver_options
.
fvectol
=
1e-6
;
options
.
solver_options
.
steptol
=
1e-14
;
options
.
system_options
.
non_ideality_tolerance
=
1e-10
;
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveSpeciation
spec_solver
(
raw_data
,
conditions
,
id_ca
,
options
);
return
spec_solver
.
get_equilibrium_curve
(
0.05
,
-
500
);
}
specmicp
::
Matrix
test_chemistry_with_al
()
{
specmicp
::
database
::
Database
thedatabase
(
"../data/cemdata_specmicp.js"
);
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
},
{
"Al[3+]"
,
"Al(OH)4[-]"
}
});
thedatabase
.
swap_components
(
swapping
);
thedatabase
.
remove_gas_phases
();
specmicp
::
RawDatabasePtr
raw_data
=
thedatabase
.
get_database
();
specmicp
::
Formulation
formulation
;
specmicp
::
scalar_t
mult
=
6.5e3
;
specmicp
::
scalar_t
m_c3s
=
mult
*
0.6
;
specmicp
::
scalar_t
m_c2s
=
mult
*
0.2
;
specmicp
::
scalar_t
m_c3a
=
mult
*
0.10
;
specmicp
::
scalar_t
m_gypsum
=
mult
*
0.10
;
specmicp
::
scalar_t
wc
=
0.8
;
specmicp
::
scalar_t
m_water
=
wc
*
1e-3
*
(
m_c3s
*
(
3
*
56.08
+
60.08
)
+
m_c2s
*
(
2
*
56.06
+
60.08
)
+
m_c3a
*
(
3
*
56.08
+
101.96
)
+
m_gypsum
*
(
56.08
+
80.06
+
2
*
18.02
)
);
formulation
.
mass_solution
=
m_water
;
formulation
.
amount_minerals
=
{
{
"C3S"
,
m_c3s
},
{
"C2S"
,
m_c2s
},
{
"C3A"
,
m_c3a
},
{
"Gypsum"
,
m_gypsum
}
};
formulation
.
minerals_to_keep
=
{
"Portlandite"
,
"CSH,jennite"
,
"CSH,tobermorite"
,
"SiO2,am"
,
"Al(OH)3,am"
,
"Monosulfoaluminate"
,
"Straetlingite"
,
"Gypsum"
,
"Ettringite"
,
};
for
(
specmicp
::
index_t
component:
raw_data
->
range_component
())
{
std
::
cout
<<
raw_data
->
labels_basis
[
component
]
<<
std
::
endl
;
}
specmicp
::
Vector
total_concentrations
=
specmicp
::
Dissolver
(
raw_data
).
dissolve
(
formulation
);
//specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp
::
index_t
id_ho
=
thedatabase
.
component_label_to_id
(
"HO[-]"
);
specmicp
::
index_t
id_ca
=
thedatabase
.
component_label_to_id
(
"Ca[2+]"
);
specmicp
::
AdimensionalSystemBC
conditions
(
total_concentrations
);
conditions
.
charge_keeper
=
id_ho
;
specmicp
::
AdimensionalSystemSolverOptions
options
;
options
.
solver_options
.
maxstep
=
20.0
;
options
.
solver_options
.
max_iter
=
100
;
options
.
solver_options
.
maxiter_maxstep
=
100
;
options
.
solver_options
.
use_crashing
=
false
;
options
.
solver_options
.
use_scaling
=
false
;
options
.
solver_options
.
factor_descent_condition
=
-
1
;
options
.
solver_options
.
factor_gradient_search_direction
=
100
;
options
.
solver_options
.
projection_min_variable
=
1e-9
;
options
.
solver_options
.
fvectol
=
1e-6
;
options
.
solver_options
.
steptol
=
1e-14
;
options
.
system_options
.
non_ideality_tolerance
=
1e-10
;
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveSpeciation
spec_solver
(
raw_data
,
conditions
,
id_ca
,
options
);
return
spec_solver
.
get_equilibrium_curve
(
0.05
,
-
500.0
);
}
void
test_eqcurve_extractor
()
{
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveExtractor
extract
(
test_chemistry_with_al
());
specmicp
::
index_t
index
=
extract
.
find_point
(
111.0
);
std
::
cout
<<
"111.0
\t
"
<<
extract
.
totsolid_concentration
(
index
)
<<
"
\t
"
<<
extract
.
totaq_concentration
(
index
)
<<
"
\t
"
<<
extract
.
porosity
(
index
)
<<
"
\t
"
<<
extract
.
diffusion_coefficient
(
index
)
<<
std
::
endl
;
}
void
test_diffeqcurve
()
{
specmicp
::
Matrix
eqcurve
=
test_chemistry_with_al
();
eqcurve
.
col
(
0
)
*=
1e-6
;
//mol/m3 -> mol/cm3
eqcurve
.
col
(
1
)
*=
1e-3
;
//mol/kg -> mol/cm3
std
::
cout
<<
eqcurve
<<
std
::
endl
;
specmicp
::
scalar_t
radius
=
3.5
;
//cm
specmicp
::
scalar_t
height
=
8.0
;
//cm
specmicp
::
scalar_t
dx
=
0.0025
;
specmicp
::
index_t
additional_nodes
=
1
;
radius
=
radius
+
additional_nodes
*
dx
;
specmicp
::
index_t
nb_nodes
=
10
+
additional_nodes
;
specmicp
::
mesh
::
Mesh1DPtr
the_mesh
=
specmicp
::
mesh
::
axisymmetric_uniform_mesh1d
(
nb_nodes
,
radius
,
dx
,
height
);
specmicp
::
dfpmsolver
::
ParabolicDriverOptions
options
;
options
.
step_tolerance
=
1e-10
;
options
.
residuals_tolerance
=
1e-8
;
options
.
quasi_newton
=
1
;
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveCoupler
solver
(
eqcurve
,
the_mesh
,
options
);
specmicp
::
scalar_t
sum_0
=
0
;
for
(
specmicp
::
index_t
node
=
0
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
{
sum_0
+=
solver
.
solid_concentrations
()(
node
)
*
the_mesh
->
get_volume_cell
(
node
);
std
::
cout
<<
the_mesh
->
get_volume_cell
(
node
)
<<
std
::
endl
;
}
specmicp
::
scalar_t
dt
=
0.5
;
specmicp
::
scalar_t
total
=
0
;
std
::
cout
<<
total
<<
"
\t
"
<<
0.0
<<
"
\t
"
<<
sum_0
<<
"
\t
"
<<
0.0
<<
std
::
endl
;
specmicp
::
index_t
k
=
0
;
while
(
total
<
0.4
)
{
solver
.
run_step
(
dt
);
total
+=
dt
/
3600
/
24
;
if
(
k
%
5000
==
0
)
{
specmicp
::
scalar_t
sum
=
0
;
for
(
specmicp
::
index_t
node
=
0
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
{
sum
+=
solver
.
solid_concentrations
()(
node
)
*
the_mesh
->
get_volume_cell
(
node
);
}
std
::
cout
<<
total
<<
"
\t
"
<<
std
::
sqrt
(
total
)
<<
"
\t
"
<<
sum
<<
"
\t
"
<<
(
sum_0
-
sum
)
/
(
1.75929e-2
)
<<
std
::
endl
;
}
++
k
;
}
//std::cout << solver.solid_concentrations() << std::endl;
}
void
test_eqcurve_solid
()
{
//specmicp::Matrix eq_curve = test_chemistry();
specmicp
::
Matrix
eq_curve
=
test_chemistry_with_al
();
eq_curve
.
col
(
0
)
*=
1e-6
;
//mol/m3 -> mol/cm3
eq_curve
.
col
(
1
)
*=
1e-3
;
//mol/kg -> mol/cm3
for
(
specmicp
::
index_t
ind
=
1
;
ind
<
eq_curve
.
rows
();
++
ind
)
{
if
(
eq_curve
(
ind
,
1
)
>=
eq_curve
(
ind
-
1
,
1
))
eq_curve
(
ind
,
1
)
=
eq_curve
(
ind
-
1
,
1
);
}
std
::
cout
<<
eq_curve
<<
std
::
endl
;
specmicp
::
scalar_t
radius
=
3.5
;
//cm
specmicp
::
scalar_t
height
=
8.0
;
//cm
specmicp
::
scalar_t
dx
=
0.005
;
specmicp
::
index_t
additional_nodes
=
1
;
radius
=
radius
+
additional_nodes
*
dx
;
specmicp
::
index_t
nb_nodes
=
50
+
additional_nodes
;
specmicp
::
mesh
::
Mesh1DPtr
the_mesh
=
specmicp
::
mesh
::
axisymmetric_uniform_mesh1d
(
nb_nodes
,
radius
,
dx
,
height
);
//specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(nb_nodes, dx, 5);
specmicp
::
dfpmsolver
::
ParabolicDriverOptions
options
;
options
.
step_tolerance
=
1e-12
;
options
.
residuals_tolerance
=
1e-6
;
options
.
sparse_solver
=
specmicp
::
SparseSolver
::
GMRES
;
//options.linesearch = specmicp::dfpmsolver::ParabolicLinesearch::Strang;
options
.
alpha
=
1.0
;
options
.
quasi_newton
=
1
;
options
.
maximum_step_length
=
10
;
specmicp
::
reactmicp
::
eqcurve
::
SolidDiffusion
program
(
the_mesh
,
eq_curve
,
{
0
,});
specmicp
::
dfpmsolver
::
ParabolicDriver
<
specmicp
::
reactmicp
::
eqcurve
::
SolidDiffusion
>
solver
(
program
);
solver
.
get_options
()
=
options
;
solver
.
set_scaling
(
specmicp
::
Vector
::
Constant
(
program
.
get_neq
(),
1e6
));
specmicp
::
Vector
variables
(
nb_nodes
);
variables
(
0
)
=
eq_curve
(
eq_curve
.
rows
()
-
10
,
0
);
for
(
specmicp
::
index_t
node
=
1
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
{
variables
(
node
)
=
eq_curve
(
10
,
0
);
}
specmicp
::
scalar_t
sum_0
=
0
;
for
(
specmicp
::
index_t
node
=
0
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
{
sum_0
+=
variables
(
node
)
*
the_mesh
->
get_volume_cell
(
node
);
std
::
cout
<<
the_mesh
->
get_volume_cell
(
node
)
<<
std
::
endl
;
}
specmicp
::
scalar_t
dt
=
10.0
;
specmicp
::
scalar_t
total
=
0
;
std
::
cout
<<
total
<<
"
\t
"
<<
0.0
<<
"
\t
"
<<
sum_0
<<
"
\t
"
<<
0.0
<<
std
::
endl
;
specmicp
::
index_t
k
=
0
;
while
(
total
<
50
)
{
//std::cout << " ==== TIMESTEP === " << std::endl;
solver
.
solve_timestep
(
dt
,
variables
);
for
(
int
node
=
0
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
if
(
variables
(
node
)
<
1e-6
)
variables
(
node
)
=
0
;
//std::cout << solver.get_perfs().nb_iterations << std::endl;
total
+=
dt
/
3600
/
24
;
if
(
k
%
5000
==
0
)
{
specmicp
::
scalar_t
sum
=
0
;
for
(
specmicp
::
index_t
node
=
0
;
node
<
the_mesh
->
nb_nodes
();
++
node
)
{
sum
+=
variables
(
node
)
*
the_mesh
->
get_volume_cell
(
node
);
}
std
::
cout
<<
total
<<
"
\t
"
<<
std
::
sqrt
(
total
)
<<
"
\t
"
<<
sum
<<
"
\t
"
<<
(
sum_0
-
sum
)
/
(
1.75929e-2
)
<<
std
::
endl
;
}
++
k
;
}
std
::
cout
<<
variables
<<
std
::
endl
;
}
void
test_interpolator
()
{
specmicp
::
Matrix
mat
(
5
,
4
);
mat
<<
1
,
1
,
1
,
1
,
2
,
1
,
2
,
0
,
3
,
1
,
3
,
-
1
,
4
,
1
,
4
,
-
2
,
5
,
1
,
5
,
-
3
;
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveExtractor
interpolator
(
mat
);
std
::
cout
<<
" "
<<
interpolator
.
slope
(
0
,
1
)
<<
" ?== "
<<
0
<<
std
::
endl
;
std
::
cout
<<
" "
<<
interpolator
.
slope
(
0
,
2
)
<<
" ?== "
<<
1
<<
std
::
endl
;
std
::
cout
<<
interpolator
.
slope
(
0
,
3
)
<<
" ?== "
<<
-
1
<<
std
::
endl
;
std
::
cout
<<
interpolator
.
find_point
(
1.5
)
<<
" ? == "
<<
0
<<
std
::
endl
;
std
::
cout
<<
interpolator
.
find_point
(
3.5
)
<<
" ? == "
<<
2
<<
std
::
endl
;
std
::
cout
<<
interpolator
.
interpolate
(
2
,
3.5
,
2
)
<<
" ? == "
<<
3.5
<<
std
::
endl
;
std
::
cout
<<
interpolator
.
interpolate
(
2
,
3.5
,
3
)
<<
" ? == "
<<
-
1.5
<<
std
::
endl
;
specmicp
::
Matrix
mat2
(
5
,
4
);
mat2
<<
5
,
1
,
1
,
1
,
4
,
1
,
2
,
0
,
3
,
1
,
3
,
-
1
,
2
,
1
,
4
,
-
2
,
1
,
1
,
5
,
-
3
;
specmicp
::
reactmicp
::
eqcurve
::
EquilibriumCurveExtractor
interpolator2
(
mat2
);
std
::
cout
<<
" "
<<
interpolator2
.
slope
(
0
,
1
)
<<
" ?== "
<<
0
<<
std
::
endl
;
std
::
cout
<<
" "
<<
interpolator2
.
slope
(
0
,
2
)
<<
" ?== "
<<
-
1
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
slope
(
0
,
3
)
<<
" ?== "
<<
+
1
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
find_point
(
1.5
)
<<
" ? == "
<<
3
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
find_point
(
3.5
)
<<
" ? == "
<<
1
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
interpolate
(
1
,
3.5
,
2
)
<<
" ? == "
<<
2.5
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
interpolate
(
1
,
3.5
,
3
)
<<
" ? == "
<<
-
0.5
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
interpolate
(
4
,
1.0
,
2
)
<<
" ? == "
<<
5
<<
std
::
endl
;
std
::
cout
<<
interpolator2
.
interpolate
(
4
,
1.0
,
3
)
<<
" ? == "
<<
-
3
<<
std
::
endl
;
}
int
main
()
{
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
//test_chemistry();
//std::cout << test_chemistry_with_al() << std::endl;
//test_eqcurve_extractor();
//test_diffeqcurve();
test_eqcurve_solid
();
test_interpolator
();
return
EXIT_SUCCESS
;
}
Event Timeline
Log In to Comment