Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77484669
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
Wed, Aug 14, 18:25
Size
19 KB
Mime Type
text/x-c
Expires
Fri, Aug 16, 18:25 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
19875905
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-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
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
;
//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
;
solver
.
get_options
().
transport_options
.
quasi_newton
=
3
;
solver
.
get_options
().
transport_options
.
linesearch
=
dfpmsolver
::
ParabolicLinesearch
::
Strang
;
// 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
();
}
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
<
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
();
}
return
EXIT_SUCCESS
;
}
int
main
()
{
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
return
leaching_neutrality
();
}
Event Timeline
Log In to Comment