Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69656372
leaching.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
Tue, Jul 2, 21:10
Size
12 KB
Mime Type
text/x-c
Expires
Thu, Jul 4, 21:10 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
18738358
Attached To
rSPECMICP SpecMiCP / ReactMiCP
leaching.cpp
View Options
#include <iostream>
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
#include "reactmicp/mesh.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include <fstream>
#include <Eigen/LU>
using
namespace
specmicp
;
using
namespace
specmicp
::
reactmicp
;
using
namespace
specmicp
::
reactmicp
::
systems
::
siasaturated
;
const
double
M_CaO
=
56.08
;
const
double
M_SiO2
=
60.09
;
const
double
M_Al2O3
=
101.96
;
const
double
M_SO3
=
80.06
;
void
compo_from_oxyde
(
Vector
&
compo_oxyde
,
Vector
&
compo_species
)
{
Eigen
::
MatrixXd
compo_in_oxyde
(
4
,
4
);
Vector
molar_mass
(
4
);
molar_mass
<<
1.0
/
M_CaO
,
1.0
/
M_SiO2
,
1.0
/
M_Al2O3
,
1.0
/
M_SO3
;
compo_oxyde
=
compo_oxyde
.
cwiseProduct
(
molar_mass
);
//C3S C2S C3A Gypsum
compo_in_oxyde
<<
3
,
2
,
3
,
1
,
// CaO
1
,
1
,
0
,
0
,
// SiO2
0
,
0
,
1
,
0
,
// Al2O3
0
,
0
,
0
,
1
;
// SO3
Eigen
::
FullPivLU
<
Eigen
::
MatrixXd
>
solver
(
compo_in_oxyde
);
compo_species
=
solver
.
solve
(
compo_oxyde
);
}
std
::
shared_ptr
<
database
::
DataContainer
>
set_database
()
{
specmicp
::
database
::
Database
database
(
"../data/cemdata_specmicp.js"
);
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
},
{
"Al[3+]"
,
"Al(OH)4[-]"
}
});
database
.
swap_components
(
swapping
);
//std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]"};
//database.keep_only_components(to_keep);
return
database
.
get_database
();
}
EquilibriumState
get_initial_state
(
std
::
shared_ptr
<
specmicp
::
database
::
DataContainer
>
database
,
Vector
&
compo
,
double
mult
)
{
std
::
shared_ptr
<
specmicp
::
ReactionPathModel
>
model
=
std
::
make_shared
<
specmicp
::
ReactionPathModel
>
();
compo
=
mult
*
compo
;
const
double
m_c3s
=
compo
(
0
);
const
double
m_c2s
=
compo
(
1
);
const
double
m_c3a
=
compo
(
2
);
const
double
m_gypsum
=
compo
(
3
);
const
double
wc
=
1.0
;
const
double
m_water
=
wc
*
(
(
3
*
M_CaO
+
M_SiO2
)
*
m_c3s
+
(
2
*
M_CaO
+
M_SiO2
)
*
m_c2s
+
(
3
*
M_CaO
+
M_Al2O3
)
*
m_c3a
+
(
M_CaO
+
M_SO3
)
*
m_gypsum
)
*
1e-3
;
model
->
amount_aqueous
=
{
{
"H2O"
,
specmicp
::
reaction_amount_t
(
m_water
/
database
->
molar_mass_basis_si
(
0
),
0.0
)},
{
"HO[-]"
,
specmicp
::
reaction_amount_t
(
0.0
,
0.0
)},
};
model
->
amount_minerals
=
{
{
"C3S"
,
specmicp
::
reaction_amount_t
(
m_c3s
,
0.0
)},
{
"C2S"
,
specmicp
::
reaction_amount_t
(
m_c2s
,
0.0
)},
{
"C3A"
,
specmicp
::
reaction_amount_t
(
m_c3a
,
0.0
)},
{
"Gypsum"
,
specmicp
::
reaction_amount_t
(
m_gypsum
,
0.0
)},
};
model
->
minerals_to_keep
=
{
"Portlandite"
,
"CSH,jennite"
,
"CSH,tobermorite"
,
"SiO2,am"
,
"Al(OH)3,am"
,
"Monosulfoaluminate"
,
"Straetlingite"
,
"Gypsum"
,
"Ettringite"
};
Eigen
::
VectorXd
x
;
specmicp
::
ReactionPathDriver
driver
(
model
,
database
,
x
);
//driver.get_options().solver_options.penalization_factor =1.0;
driver
.
get_options
().
ncp_function
=
specmicp
::
micpsolver
::
NCPfunction
::
penalizedFB
;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver
.
get_options
().
solver_options
.
use_scaling
=
false
;
driver
.
get_options
().
solver_options
.
max_factorization_step
=
4
;
driver
.
get_options
().
solver_options
.
factor_gradient_search_direction
=
200
;
driver
.
get_options
().
solver_options
.
maxstep
=
50
;
driver
.
get_options
().
solver_options
.
maxiter_maxstep
=
50
;
driver
.
get_options
().
solver_options
.
max_iter
=
100
;
driver
.
get_options
().
allow_restart
=
true
;
specmicp
::
micpsolver
::
MiCPPerformance
perf
=
driver
.
one_step
(
x
);
if
(
perf
.
return_code
!=
specmicp
::
micpsolver
::
MiCPSolverReturnCode
::
ResidualMinimized
)
throw
std
::
runtime_error
(
"Error - system cannot be solved"
);
return
driver
.
get_current_solution
();
}
EquilibriumState
get_blank_state
(
std
::
shared_ptr
<
specmicp
::
database
::
DataContainer
>
database
)
{
std
::
shared_ptr
<
specmicp
::
ReactionPathModel
>
model
=
std
::
make_shared
<
specmicp
::
ReactionPathModel
>
();
const
double
m_c3s
=
1e-7
;
const
double
m_c2s
=
1e-7
;
const
double
m_c3a
=
1e-7
;
const
double
m_gypsum
=
1e-7
;
const
double
m_water
=
1
;
model
->
amount_aqueous
=
{
{
"H2O"
,
specmicp
::
reaction_amount_t
(
m_water
/
database
->
molar_mass_basis_si
(
0
),
0.0
)},
{
"HO[-]"
,
specmicp
::
reaction_amount_t
(
-
0.0
,
0.0
)},
{
"SiO(OH)3[-]"
,
specmicp
::
reaction_amount_t
(
0.0
,
0.0
)}
};
model
->
amount_minerals
=
{
{
"C3S"
,
specmicp
::
reaction_amount_t
(
m_c3s
,
0.0
)},
{
"C2S"
,
specmicp
::
reaction_amount_t
(
m_c2s
,
0.0
)},
{
"C3A"
,
specmicp
::
reaction_amount_t
(
m_c3a
,
0.0
)},
{
"Gypsum"
,
specmicp
::
reaction_amount_t
(
m_gypsum
,
0.0
)},
{
"SiO2,am"
,
specmicp
::
reaction_amount_t
(
0.0001
,
0.0
)}
};
Eigen
::
VectorXd
x
;
specmicp
::
ReactionPathDriver
driver
(
model
,
database
,
x
);
//driver.get_options().solver_options.penalization_factor =1.0;
driver
.
get_options
().
ncp_function
=
specmicp
::
micpsolver
::
NCPfunction
::
penalizedFB
;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver
.
get_options
().
solver_options
.
use_scaling
=
false
;
driver
.
get_options
().
solver_options
.
max_factorization_step
=
1
;
driver
.
get_options
().
solver_options
.
factor_gradient_search_direction
=
50
;
driver
.
get_options
().
solver_options
.
maxstep
=
50
;
driver
.
get_options
().
solver_options
.
maxiter_maxstep
=
50
;
driver
.
get_options
().
solver_options
.
max_iter
=
100
;
driver
.
get_options
().
allow_restart
=
true
;
specmicp
::
micpsolver
::
MiCPPerformance
perf
=
driver
.
one_step
(
x
);
if
(
perf
.
return_code
!=
specmicp
::
micpsolver
::
MiCPSolverReturnCode
::
ResidualMinimized
)
throw
std
::
runtime_error
(
"Error - blank system cannot be solved"
);
return
driver
.
get_current_solution
();
}
double
get_dt
(
std
::
shared_ptr
<
SaturatedDiffusionTransportParameters
>
parameters
,
std
::
shared_ptr
<
mesh
::
Mesh1D
>
themesh
)
{
//const scalar_t max_d = parameters->diffusion_coefficient(1);
//scalar_t dt = std::pow(themesh->get_dx(0), 2)/max_d;
//if (dt<3600*24) dt=3600*24;
//if (dt>2*3600*24) dt=2*3600*24;
return
3600
*
24
/
5
;
}
int
main
()
{
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
Vector
compo_oxyde
(
4
),
compo_species
(
4
);
compo_oxyde
<<
62.9
,
20.6
,
5.8
,
3.1
;
compo_from_oxyde
(
compo_oxyde
,
compo_species
);
std
::
shared_ptr
<
database
::
DataContainer
>
thedatabase
=
set_database
();
scalar_t
radius
=
3.5
;
//cm
scalar_t
height
=
8
;
//cm
index_t
nb_nodes
=
36
;
mesh
::
Mesh1DPtr
themesh
=
mesh
::
axisymmetric_uniform_mesh1d
(
nb_nodes
,
radius
,
height
);
std
::
shared_ptr
<
SaturatedDiffusionTransportParameters
>
parameters
=
std
::
make_shared
<
SaturatedDiffusionTransportParameters
>
(
nb_nodes
,
1e4
*
std
::
exp
(
9.95
*
0.25
-
29.08
),
0.25
);
parameters
->
diffusion_coefficient
(
0
)
=
2.2e-5
;
parameters
->
porosity
(
0
)
=
1.0
;
EquilibriumState
initial_state
=
get_initial_state
(
thedatabase
,
compo_species
,
0.5
);
SIABoundaryConditions
bcs
(
nb_nodes
);
bcs
.
list_initial_states
.
push_back
(
initial_state
);
bcs
.
list_initial_states
.
push_back
(
get_blank_state
(
thedatabase
));
bcs
.
bs_types
[
0
]
=
BoundaryConditionType
::
FixedComposition
;
bcs
.
initial_states
[
0
]
=
1
;
database
::
Database
dbhandler
(
thedatabase
);
index_t
id_ca
=
dbhandler
.
component_label_to_id
(
"Ca[2+]"
);
index_t
id_si
=
dbhandler
.
component_label_to_id
(
"SiO(OH)3[-]"
);
std
::
ofstream
output_ca
,
output_aqca
,
output_si
,
output_poro
;
output_ca
.
open
(
"out_ca.dat"
);
output_aqca
.
open
(
"out_aqca.dat"
);
output_si
.
open
(
"out_si.dat"
);
output_poro
.
open
(
"out_poro.ca"
);
std
::
vector
<
std
::
ofstream
>
mineral_out
(
thedatabase
->
nb_mineral
);
for
(
int
mineral:
thedatabase
->
range_mineral
())
{
mineral_out
[
mineral
].
open
(
"out_mineral_"
+
std
::
to_string
(
mineral
)
+
".dat"
);
}
SIASaturatedReactiveTransportSolver
solver
(
themesh
,
thedatabase
,
parameters
);
solver
.
apply_boundary_conditions
(
bcs
);
solver
.
get_variables
().
scale_volume
();
solver
.
use_sia
(
100
);
solver
.
get_options
().
residual_tolerance
=
1e-3
;
solver
.
get_options
().
threshold_update_disable_speciation
=
1e-16
;
solver
.
get_options
().
threshold_update_transport
=
1e-10
;
// Transport
solver
.
get_options
().
transport_options
.
residuals_tolerance
=
1e-5
;
// Speciation
solver
.
get_options
().
speciation_options
.
solver_options
.
use_crashing
=
false
;
solver
.
get_options
().
speciation_options
.
solver_options
.
use_scaling
=
false
;
solver
.
get_options
().
speciation_options
.
solver_options
.
max_iter
=
100
;
solver
.
get_options
().
speciation_options
.
solver_options
.
maxstep
=
50
;
std
::
cout
<<
"id SiO2,am : "
<<
dbhandler
.
mineral_label_to_id
(
"SiO2,am"
)
<<
std
::
endl
;
//thedatabase->stability_mineral(
// dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve;
//double dt = 5e4;
double
total
=
0
;
uindex_t
k
=
1
;
while
(
total
<
365
*
24
*
3600
)
{
double
dt
=
get_dt
(
parameters
,
themesh
);
total
+=
dt
;
std
::
cout
<<
" iterations : "
<<
k
<<
" - dt : "
<<
dt
<<
std
::
endl
;
SIASaturatedReactiveTransportSolverReturnCode
retcode
=
solver
.
solve_timestep
(
dt
);
if
(
retcode
!=
SIASaturatedReactiveTransportSolverReturnCode
::
Success
)
throw
std
::
runtime_error
(
"Failed to converge : "
+
std
::
to_string
((
int
)
retcode
));
SIASaturatedVariables
&
variables
=
solver
.
get_variables
();
//std::cout << variables.total_mobile_concentrations() << std::endl;
// compute porosity
for
(
index_t
node:
themesh
->
range_nodes
())
{
double
volume
=
variables
.
volume_minerals
(
node
);
double
cell_volume
=
themesh
->
get_volume_cell
(
node
);
specmicp_assert
(
volume
<
cell_volume
);
double
porosity
=
(
cell_volume
-
volume
*
1e6
)
/
cell_volume
;
parameters
->
porosity
(
node
)
=
porosity
;
parameters
->
diffusion_coefficient
(
node
)
=
porosity
>
0.6
?
1e4
*
std
::
exp
(
9.95
*
0.6
-
29.08
)
:
1e4
*
std
::
exp
(
9.95
*
porosity
-
29.08
);
//std::cout << parameters->diffusion_coefficient(node);
}
parameters
->
diffusion_coefficient
(
0
)
=
parameters
->
diffusion_coefficient
(
1
);
//std::cout << variables.speciation_variables() << std::endl;
if
(
k
<
10
or
k
%
5
==
0
)
{
output_ca
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
output_aqca
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
output_si
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
output_poro
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
for
(
index_t
node:
themesh
->
range_nodes
())
{
output_ca
<<
" "
<<
variables
.
nodal_component_update_total_amount
(
node
,
id_ca
);
output_aqca
<<
" "
<<
variables
.
total_mobile_concentration
(
node
,
id_ca
);
output_si
<<
" "
<<
variables
.
nodal_component_update_total_amount
(
node
,
id_si
);
double
volume
=
variables
.
volume_minerals
(
node
);
double
cell_volume
=
themesh
->
get_volume_cell
(
node
);
output_poro
<<
" "
<<
(
cell_volume
-
volume
*
1e6
)
/
cell_volume
;
}
output_ca
<<
std
::
endl
;
output_aqca
<<
std
::
endl
;
output_si
<<
std
::
endl
;
output_poro
<<
std
::
endl
;
for
(
index_t
mineral:
thedatabase
->
range_mineral
())
{
mineral_out
[
mineral
]
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
for
(
index_t
node:
themesh
->
range_nodes
())
{
mineral_out
[
mineral
]
<<
" "
<<
variables
.
mineral_amount
(
node
,
mineral
);
}
mineral_out
[
mineral
]
<<
std
::
endl
;
}
}
++
k
;
}
output_ca
.
close
();
output_aqca
.
close
();
output_si
.
close
();
output_poro
.
close
();
for
(
index_t
mineral:
thedatabase
->
range_mineral
())
{
mineral_out
[
mineral
].
close
();
}
}
Event Timeline
Log In to Comment