Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78697325
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
Thu, Aug 22, 11:28
Size
23 KB
Mime Type
text/x-c
Expires
Sat, Aug 24, 11:28 (2 d)
Engine
blob
Format
Raw Data
Handle
20087602
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 "reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.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-8;
//const double m_c2s = 1e-8;
//const double m_c3a = 10e-8;
//const double m_gypsum = 10e-8;
//const double m_water = 1.0;
//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(database->nb_component+database->nb_mineral);
//x.setZero();
//x(0) = 1.0;
//x.segment(1, database->nb_component-1).setConstant(-6);
//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");
//
//EquilibriumState solution = driver.get_current_solution();
//std::cout << solution.main_variables() << std::endl;
//solution.remove_minerals();
Eigen
::
VectorXd
variables
(
database
->
nb_component
+
database
->
nb_mineral
);
variables
.
setZero
();
variables
(
0
)
=
1.0
;
for
(
index_t
component
=
1
;
component
<
database
->
nb_component
;
++
component
)
variables
(
component
)
=
-
11
;
scalar_t
ionic_strength
=
0.0
;
Eigen
::
VectorXd
secondary
(
database
->
nb_aqueous
);
Eigen
::
VectorXd
loggamma
(
database
->
nb_component
+
database
->
nb_aqueous
);
secondary
.
setZero
();
loggamma
.
setZero
();
EquilibriumState
solution
(
variables
,
secondary
,
loggamma
,
ionic_strength
,
database
);
return
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
/
3
;
if
(
dt
>
5.0
)
return
5.0
;
return
5.0
;
//return 10.0;
//return dt;
//if (dt<3600*24) dt=3600*24;
//if (dt>2*3600*24) dt=2*3600*24;
//return 900;
}
int
leaching
()
{
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.0
;
//cm
scalar_t
dx
=
0.005
;
index_t
additional_nodes
=
1
;
radius
=
radius
+
additional_nodes
*
dx
;
index_t
nb_nodes
=
50
;
mesh
::
Mesh1DPtr
themesh
=
mesh
::
axisymmetric_uniform_mesh1d
(
nb_nodes
,
radius
,
dx
,
height
);
std
::
shared_ptr
<
SaturatedDiffusionTransportParameters
>
parameters
=
std
::
make_shared
<
SaturatedDiffusionTransportParameters
>
(
nb_nodes
,
1e4
*
std
::
exp
(
9.95
*
0.15
-
29.08
),
0.15
);
parameters
->
diffusion_coefficient
(
0
)
=
2.2e-5
;
// cm^2/s
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;
for
(
index_t
node
=
0
;
node
<
additional_nodes
;
++
node
)
bcs
.
initial_states
[
node
]
=
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[-]"
);
index_t
id_oh
=
dbhandler
.
component_label_to_id
(
"HO[-]"
);
for
(
index_t
mineral:
thedatabase
->
range_mineral
())
{
std
::
cout
<<
thedatabase
->
labels_minerals
[
mineral
]
<<
std
::
endl
;
}
std
::
cout
<<
"Volume first cell "
<<
themesh
->
get_volume_cell
(
additional_nodes
)
<<
std
::
endl
;
//return EXIT_SUCCESS;
std
::
ofstream
output_ca
,
output_aqca
,
output_si
,
output_poro
,
output_ph
;
std
::
ofstream
output_saca
,
output_sasi
;
output_ca
.
open
(
"out_ca.dat"
);
output_aqca
.
open
(
"out_aqca.dat"
);
output_si
.
open
(
"out_si.dat"
);
output_ph
.
open
(
"out_ph.dat"
);
output_poro
.
open
(
"out_poro.dat"
);
output_saca
.
open
(
"out_sa_ca.dat"
);
output_sasi
.
open
(
"out_sa_si.dat"
);
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_snia();
solver
.
use_sia
(
100
);
solver
.
get_options
().
residual_tolerance
=
1e-4
;
solver
.
get_options
().
threshold_update_disable_speciation
=
1e-16
;
solver
.
get_options
().
threshold_update_transport
=
1e-10
;
solver
.
get_options
().
good_enough_transport_residual_tolerance
=
1e-3
;
solver
.
get_options
().
relative_update
=
1e-4
;
solver
.
get_options
().
use_consistent_velocity
=
true
;
// Transport
solver
.
get_options
().
transport_options
.
residuals_tolerance
=
1e-6
;
solver
.
get_options
().
transport_options
.
maximum_iterations
=
500
;
solver
.
get_options
().
transport_options
.
step_tolerance
=
1e-8
;
solver
.
get_options
().
transport_options
.
threshold_stationary_point
=
1e-5
;
solver
.
get_options
().
transport_options
.
quasi_newton
=
5
;
solver
.
get_options
().
transport_options
.
sparse_solver
=
SparseSolver
::
GMRES
;
//solver.get_options().transport_options.alpha = 0.0;
solver
.
get_options
().
transport_options
.
linesearch
=
dfpmsolver
::
ParabolicLinesearch
::
Strang
;
// Speciation
solver
.
get_options
().
speciation_options
.
solver_options
.
fvectol
=
1e-10
;
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
=
50
;
solver
.
get_options
().
speciation_options
.
solver_options
.
maxstep
=
50
;
solver
.
get_options
().
speciation_options
.
solver_options
.
non_monotone_linesearch
=
true
;
solver
.
get_options
().
speciation_options
.
solver_options
.
threshold_cycling_linesearch
=
1e-8
;
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;
SIASaturatedVariables
&
first_variables
=
solver
.
get_variables
();
std
::
cout
<<
first_variables
.
total_mobile_concentrations
()
<<
std
::
endl
;
//double dt = 5e4;
double
total
=
0
;
uindex_t
k
=
0
;
scalar_t
dt
=
50.0
;
while
(
total
<
365
*
24
*
3600
)
{
total
+=
dt
;
std
::
cout
<<
" iterations : "
<<
k
<<
" - dt : "
<<
dt
<<
std
::
endl
;
SIASaturatedReactiveTransportSolverReturnCode
retcode
=
solver
.
solve_timestep
(
dt
);
std
::
cout
<<
" Return code : "
<<
(
int
)
retcode
<<
" - Nb iter : "
<<
solver
.
get_perfs
().
nb_iterations
<<
std
::
endl
<<
" - Nb iter transport : "
<<
solver
.
get_perfs
().
nb_iterations_transport
<<
std
::
endl
<<
" - Nb iter chemistry : "
<<
solver
.
get_perfs
().
nb_iterations_chemistry
<<
std
::
endl
;
if
(
retcode
!=
SIASaturatedReactiveTransportSolverReturnCode
::
Success
)
{
std
::
cout
<<
"Current residual : "
<<
solver
.
get_perfs
().
current_residual
<<
std
::
endl
;
throw
std
::
runtime_error
(
"Failed to converge : "
+
std
::
to_string
((
int
)
retcode
));
//dt = dt/1.5;
//goto start_again;
}
//if (solver.get_perfs().nb_iterations <= 10) dt*=1.1;
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.92
?
2.219e-5
:
1e4
*
std
::
exp
(
9.95
*
porosity
-
29.08
);
}
std
::
cout
<<
variables
.
component_concentration
(
1
,
id_ca
)
<<
std
::
endl
;
parameters
->
diffusion_coefficient
(
0
)
=
2.219e-5
;
// parameters->diffusion_coefficient(1);
//std::cout << variables.speciation_variables() << std::endl;
if
(
k
==
1
or
(
k
<
1000
and
k
%
100
==
0
)
or
k
%
1000
==
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
);
output_ph
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
output_saca
<<
total
<<
" "
<<
total
/
(
3600.0
*
24.0
);
output_sasi
<<
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
);
output_ph
<<
" "
<<
14
+
variables
.
log_gamma_component
(
node
,
id_oh
)
+
variables
.
log_component_concentration
(
node
,
id_oh
);
double
volume
=
variables
.
volume_minerals
(
node
);
double
cell_volume
=
themesh
->
get_volume_cell
(
node
);
output_poro
<<
" "
<<
(
cell_volume
-
volume
*
1e6
)
/
cell_volume
;
output_saca
<<
" "
<<
variables
.
component_concentration_in_minerals
(
node
,
id_ca
);
output_sasi
<<
" "
<<
variables
.
component_concentration_in_minerals
(
node
,
id_si
);
}
output_ca
<<
std
::
endl
;
output_aqca
<<
std
::
endl
;
output_si
<<
std
::
endl
;
output_poro
<<
std
::
endl
;
output_ph
<<
std
::
endl
;
output_saca
<<
std
::
endl
;
output_sasi
<<
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
();
output_saca
.
close
();
output_sasi
.
close
();
for
(
index_t
mineral:
thedatabase
->
range_mineral
())
{
mineral_out
[
mineral
].
close
();
}
return
EXIT_SUCCESS
;
}
int
leaching_neutrality
()
{
std
::
shared_ptr
<
database
::
DataContainer
>
thedatabase
=
set_database
();
Vector
compo_oxyde
(
4
),
compo_species
(
4
);
compo_oxyde
<<
62.9
,
20.6
,
5.8
,
3.1
;
compo_from_oxyde
(
compo_oxyde
,
compo_species
);
scalar_t
radius
=
3.5
;
//cm
scalar_t
height
=
8
;
//cm
index_t
nb_nodes
=
36
;
database
::
Database
dbhandler
(
thedatabase
);
//mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height);
mesh
::
Mesh1DPtr
themesh
=
mesh
::
uniform_mesh1d
(
nb_nodes
,
0.1
,
5
);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std
::
shared_ptr
<
SaturatedNeutralityDiffusionTransportParameters
>
parameters
=
std
::
make_shared
<
SaturatedNeutralityDiffusionTransportParameters
>
(
nb_nodes
,
// nb_nodes,
thedatabase
->
nb_component
,
// nb_component,
thedatabase
->
nb_aqueous
,
// nb_aqueous,
2e-5
,
// the_diffusion_coefficient,
0.25
,
// porosity
0.5e-3
//reduction_factor
);
//parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5;
//parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5;
parameters
->
diff_coeff_secondary
(
dbhandler
.
aqueous_label_to_id
(
"Si(OH)4"
))
=
1e-7
;
parameters
->
diff_coeff_component
(
dbhandler
.
component_label_to_id
(
"SiO(OH)3[-]"
))
=
1e-7
;
parameters
->
porosity
(
0
)
=
1.0
;
parameters
->
reduction_factor
(
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
;
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"
);
}
SIASaturatedReactiveTransportNeutralitySolver
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-6
;
// Transport
solver
.
get_options
().
transport_options
.
residuals_tolerance
=
1e-5
;
solver
.
get_options
().
transport_options
.
maximum_iterations
=
100
;
solver
.
get_options
().
transport_options
.
quasi_newton
=
2
;
solver
.
get_options
().
transport_options
.
maximum_step_length
=
1000
;
solver
.
get_options
().
transport_options
.
linesearch
=
dfpmsolver
::
ParabolicLinesearch
::
Strang
;
solver
.
get_options
().
transport_options
.
max_iterations_at_max_length
=
100
;
// 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
;
std
::
cout
<<
"Basis :"
;
for
(
index_t
component:
thedatabase
->
range_component
())
{
std
::
cout
<<
" "
<<
thedatabase
->
labels_basis
[
component
];
}
std
::
cout
<<
std
::
endl
;
std
::
cout
<<
"Solid :"
;
for
(
index_t
mineral:
thedatabase
->
range_mineral
())
{
std
::
cout
<<
" "
<<
thedatabase
->
labels_minerals
[
mineral
];
}
std
::
cout
<<
std
::
endl
;
while
(
total
<
10
*
24
*
3600
)
{
double
dt
=
3600
/
2
;
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
->
reduction_factor
(
node
)
=
porosity
>
0.6
?
1e4
*
std
::
exp
(
9.95
*
0.6
-
29.08
)
/
2e-5
:
1e4
*
std
::
exp
(
9.95
*
porosity
-
29.08
)
/
2e-5
;
//std::cout << parameters->diffusion_coefficient(node);
}
parameters
->
reduction_factor
(
0
)
=
parameters
->
reduction_factor
(
1
);
//std::cout << variables.speciation_variables() << std::endl;
if
(
k
<
2
or
k
%
100
==
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
();
}
return
EXIT_SUCCESS
;
}
int
main
()
{
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
return
leaching
();
//return leaching_neutrality();
}
Event Timeline
Log In to Comment