Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81222400
carbonation.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, 20:21
Size
26 KB
Mime Type
text/x-c++
Expires
Fri, Sep 6, 20:21 (2 d)
Engine
blob
Format
Raw Data
Handle
20530483
Attached To
rSPECMICP SpecMiCP / ReactMiCP
carbonation.cpp
View Options
// ============================================
// Saturated carbonation of a cement past
// =============================================
// This file can be used to learn how to use ReactMiCP
// Include ReactMiCP
#include "reactmicp.hpp"
// Other includes that would be needed later :
#include <Eigen/QR>
#include "physics/laws.hpp"
#include "utils/timer.hpp"
#include <functional>
#include <iostream>
#include <fstream>
// we use the following namespaces
using
namespace
specmicp
;
using
namespace
reactmicp
;
using
namespace
reactmicp
::
systems
;
using
namespace
specmicp
::
reactmicp
::
systems
::
satdiff
;
using
namespace
reactmicp
::
solver
;
// We declare some functions that we will be using :
// Database
// ========
// Parse and Prepare the database
RawDatabasePtr
get_database
(
const
std
::
string
&
path_to_database
);
// Boundary and initial conditions
// ===============================
// Return the chemistry constraints for the inital condition
AdimensionalSystemConstraints
get_cement_initial_constraints
(
RawDatabasePtr
the_database
,
const
units
::
UnitsSet
&
units_set
);
// return the chemistry constraints for the boundary condition
AdimensionalSystemConstraints
get_bc_initial_constraints
(
RawDatabasePtr
the_database
,
const
units
::
UnitsSet
&
units_set
);
Vector
initial_compo
(
const
Vector
&
publi
);
void
compo_from_oxyde
(
Vector
&
compo_oxyde
,
Vector
&
compo_species
);
// The mesh
// ========
// There is two version
// a uniform mesh
mesh
::
Mesh1DPtr
get_mesh
(
scalar_t
dx
);
// and a non uniform mesh
mesh
::
Mesh1DPtr
get_mesh
();
// Options
// =======
// Speciation solver options
AdimensionalSystemSolverOptions
get_specmicp_options
();
// Transport options
void
set_transport_options
(
dfpmsolver
::
ParabolicDriverOptions
&
transport_options
);
// Reactive Solver options
void
set_reactive_solver_options
(
reactmicp
::
solver
::
ReactiveTransportOptions
&
reactive_options
);
// Upscaling
// =========
// Upscaling is the stagger in charge of finding the porosity and transport properties
// This is a class that is dependant on the problem to solve and should be defined by the user
// The following class is an example on how to define such a class
// The Upscaling class should inherit from UpscalingStaggerBase
class
CarboUspcaling
:
public
solver
::
UpscalingStaggerBase
{
public
:
// Initiation
// This is not call automatically so can be adjusted by the user
CarboUspcaling
(
RawDatabasePtr
the_database
,
units
::
UnitsSet
the_units
)
:
m_data
(
the_database
),
m_units_set
(
the_units
),
m_dt
(
HUGE_VAL
)
{
}
// Initialize the stagger at the beginning of the computation
// The porosity and transport properties should be initialized here
void
initialize
(
VariablesBasePtr
var
)
{
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
var
);
true_var
->
upscaling_variables
().
setZero
();
for
(
index_t
node
=
0
;
node
<
true_var
->
nb_nodes
();
++
node
)
{
upscaling_one_node
(
node
,
true_var
);
}
}
// Initialize the stagger at the beginning of a timestep
// Typically do nothing but erase the "dot" variable (velocity such as the vel_porosity)
void
initialize_timestep
(
scalar_t
dt
,
VariablesBasePtr
var
)
{
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
var
);
m_dt
=
dt
;
for
(
index_t
node
=
1
;
node
<
true_var
->
nb_nodes
();
++
node
)
{
upscaling_one_node
(
node
,
true_var
);
true_var
->
vel_porosity
(
node
)
=
0.0
;
}
}
//! This is the main function called during a timestep
StaggerReturnCode
restart_timestep
(
VariablesBasePtr
var
)
{
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
var
);
for
(
index_t
node
=
1
;
node
<
true_var
->
nb_nodes
();
++
node
)
{
upscaling_one_node
(
node
,
true_var
);
}
return
StaggerReturnCode
::
ResidualMinimized
;
// This is the value that should be returned if everything is ok
}
// A function that return the diffusion coefficient when porosity is equal to 'porosity'
scalar_t
coeff_diffusion_law
(
scalar_t
porosity
)
{
static
constexpr
scalar_t
de_0
=
2.726e-12
;
//1e-11;
static
constexpr
scalar_t
factor
=
1e4
;
static
constexpr
scalar_t
porosity_r
=
0.02
;
const
scalar_t
ratio
=
((
porosity
-
porosity_r
)
/
(
0.25
-
porosity_r
));
const
scalar_t
d_e
=
factor
*
de_0
*
std
::
pow
(
ratio
,
3.32
);
return
std
::
min
(
d_e
,
1e4
*
1e-10
);
}
// Compute the upscaling for 'node'
void
upscaling_one_node
(
index_t
node
,
SaturatedVariablesPtr
true_var
)
{
// AdimensionalSystemSolutionExtractor is the class to use to
// extract information from a SpecMiCP solution
// To obtain correct information the correct units must be used
AdimensionalSystemSolutionExtractor
extractor
(
true_var
->
equilibrium_solution
(
node
),
m_data
,
m_units_set
);
// We can obtain the porosity very easily :
scalar_t
porosity
=
extractor
.
porosity
();
true_var
->
vel_porosity
(
node
)
+=
(
porosity
-
true_var
->
porosity
(
node
))
/
m_dt
;
true_var
->
porosity
(
node
)
=
porosity
;
true_var
->
diffusion_coefficient
(
node
)
=
coeff_diffusion_law
(
porosity
);
}
private
:
RawDatabasePtr
m_data
;
units
::
UnitsSet
m_units_set
;
index_t
m_id_cshj
;
scalar_t
m_initial_sat_cshj
;
scalar_t
m_dt
;
};
// ==============
// Output
// ==============
// Many information is generated throughout a computation and we can't save everything
// The reactmicp runner can call a function to output the desired information
// Here the output function is a member function for convenience
class
OutputPolicy
{
public
:
// Initialise the output
// SimulationInformation contain basic information from the simulation
OutputPolicy
(
const
SimulationInformation
&
simul_info
,
mesh
::
Mesh1DPtr
the_mesh
,
RawDatabasePtr
the_database
,
units
::
UnitsSet
the_units
)
:
m_mesh
(
the_mesh
),
m_database
(
the_database
),
m_units
(
the_units
)
{
// Open the files
// Save the pH
out_c_ph
.
open
(
simul_info
.
complete_filepath
(
"ph"
,
suffix
));
// Save the total aqueous concentration of HCO3[-]
out_c_hco3
.
open
(
simul_info
.
complete_filepath
(
"c_hco3"
,
suffix
));
// save the total solid concentration of s_ca
out_s_ca
.
open
(
simul_info
.
complete_filepath
(
"s_ca"
,
suffix
));
// Save the volume fraction of calcite
out_calcite
.
open
(
simul_info
.
complete_filepath
(
"calcite"
,
suffix
));
// Save the porosity
out_poro
.
open
(
simul_info
.
complete_filepath
(
"poro"
,
suffix
));
// save useful indices from the database
m_id_hco3
=
m_database
->
get_id_component
(
"HCO3[-]"
);
m_id_ca
=
m_database
->
get_id_component
(
"Ca[2+]"
);
m_id_calcite
=
m_database
->
get_id_mineral
(
"Calcite"
);
}
void
output
(
scalar_t
time
,
VariablesBasePtr
variables
)
{
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
variables
);
out_c_ph
<<
time
;
out_c_hco3
<<
time
;
out_s_ca
<<
time
;
out_calcite
<<
time
;
out_poro
<<
time
;
for
(
index_t
node:
m_mesh
->
range_nodes
())
{
// Again AdimensionalSystemSolutionExtractor is used to obtain information about the speciation system
AdimensionalSystemSolutionExtractor
extractor
(
true_var
->
equilibrium_solution
(
node
),
m_database
,
m_units
);
out_c_ph
<<
"
\t
"
<<
extractor
.
pH
();
out_c_hco3
<<
"
\t
"
<<
true_var
->
aqueous_concentration
(
node
,
m_id_hco3
);
out_s_ca
<<
"
\t
"
<<
true_var
->
solid_concentration
(
node
,
m_id_ca
);
out_calcite
<<
"
\t
"
<<
extractor
.
volume_fraction_mineral
(
m_id_calcite
);
out_poro
<<
"
\t
"
<<
extractor
.
porosity
();
}
out_c_hco3
<<
std
::
endl
;
out_c_ph
<<
std
::
endl
;
out_s_ca
<<
std
::
endl
;
out_calcite
<<
std
::
endl
;
out_poro
<<
std
::
endl
;
}
private
:
mesh
::
Mesh1DPtr
m_mesh
;
RawDatabasePtr
m_database
;
units
::
UnitsSet
m_units
;
std
::
string
suffix
{
"dat"
};
std
::
ofstream
out_c_ph
;
std
::
ofstream
out_c_hco3
;
std
::
ofstream
out_s_ca
;
std
::
ofstream
out_calcite
;
std
::
ofstream
out_poro
;
index_t
m_id_hco3
;
index_t
m_id_ca
;
index_t
m_id_calcite
;
};
// ===================
// MAIN
// ===================
int
main
()
{
// initialize eigen
// This is needed if OpenMP is used
Eigen
::
initParallel
();
// Setup the log
// The output will be on cerr
specmicp
::
logger
::
ErrFile
::
stream
()
=
&
std
::
cerr
;
// And only Warning and error messages will be printed
specmicp
::
stdlog
::
ReportLevel
()
=
specmicp
::
logger
::
Warning
;
// The database
// -------------
// Obtain the database
RawDatabasePtr
the_database
=
get_database
(
"../data/cemdata_specmicp.js"
);
// SpecMiCP options
// ----------------
// Set the options
AdimensionalSystemSolverOptions
specmicp_options
=
get_specmicp_options
();
units
::
UnitsSet
units_set
=
specmicp_options
.
units_set
;
// The mesh
// --------
// Obtain the mesh
mesh
::
Mesh1DPtr
the_mesh
=
get_mesh
(
0.0050
);
// uniform
//mesh::Mesh1DPtr the_mesh = get_mesh(); // non uniform
index_t
nb_nodes
=
the_mesh
->
nb_nodes
();
// just boilerplate to make things easier
// Print the mesh
io
::
print_mesh
(
"out_carbo_mesh.dat"
,
the_mesh
,
units
::
LengthUnit
::
centimeter
);
// Initial and boundary conditions
// --------------------------------
// Obtain the initial constraints
// Note : the database is modified at this point
// Only the relevant components are kept at this point
AdimensionalSystemConstraints
cement_constraints
=
get_cement_initial_constraints
(
the_database
,
units_set
);
// Freeze the database
the_database
->
freeze_db
();
// Obtain the initial condition
AdimensionalSystemSolution
cement_solution
=
solve_equilibrium
(
the_database
,
cement_constraints
,
specmicp_options
);
if
(
not
the_database
->
is_valid
())
throw
std
::
runtime_error
(
"The database has been changed during the initial condition computation"
);
// Just to check that everything is ok !
AdimensionalSystemSolutionExtractor
extractor
(
cement_solution
,
the_database
,
units_set
);
std
::
cout
<<
"Porosity : "
<<
extractor
.
porosity
()
<<
std
::
endl
;
std
::
cout
<<
"Water volume fraction "
<<
extractor
.
volume_fraction_water
()
<<
std
::
endl
;
//return EXIT_SUCCESS;
// Obtain the boundary conditions
AdimensionalSystemConstraints
bc_constraints
=
get_bc_initial_constraints
(
the_database
,
specmicp_options
.
units_set
);
AdimensionalSystemSolution
bc_solution
=
solve_equilibrium
(
the_database
,
bc_constraints
,
specmicp_options
);
if
(
not
the_database
->
is_valid
())
throw
std
::
runtime_error
(
"The database has been changed during the boundary conditions computation"
);
// The constraints
// --------------
// Boundary condition
// 'list_fixed_nodes' lists the node with fixed composition
std
::
vector
<
index_t
>
list_fixed_nodes
=
{
0
};
// list_constraints lists the different speciation constraints in the system
// Note : the total concentrations are computed from the initial solution,
// and thus different total concentrations do not constitute a reason to create different constraints
std
::
vector
<
specmicp
::
AdimensionalSystemConstraints
>
list_constraints
=
{
cement_constraints
,
bc_constraints
};
// index_constraints links a node to the set of constraints that will be used
std
::
vector
<
int
>
index_constraints
(
nb_nodes
,
0
);
index_constraints
[
0
]
=
1
;
// Boundary conditions
// This is the list of initial composition
std
::
vector
<
specmicp
::
AdimensionalSystemSolution
>
list_initial_composition
=
{
cement_solution
,
bc_solution
};
// This list links a node to its initial conditions
std
::
vector
<
int
>
index_initial_state
(
nb_nodes
,
0
);
index_initial_state
[
0
]
=
1
;
// boundary condition
// Variables
// -----------
// Create the main set of variables
// They will be initialized using the list of initial composition
satdiff
::
SaturatedVariablesPtr
variables
=
satdiff
::
init_variables
(
the_mesh
,
the_database
,
units_set
,
list_fixed_nodes
,
list_initial_composition
,
index_initial_state
);
// Staggers
// ---------
// This section initialize the staggers
std
::
shared_ptr
<
EquilibriumStagger
>
equilibrium_stagger
=
std
::
make_shared
<
satdiff
::
EquilibriumStagger
>
(
list_constraints
,
index_constraints
,
specmicp_options
);
std
::
shared_ptr
<
SaturatedTransportStagger
>
transport_stagger
=
std
::
make_shared
<
satdiff
::
SaturatedTransportStagger
>
(
variables
,
list_fixed_nodes
);
set_transport_options
(
transport_stagger
->
options_solver
());
UpscalingStaggerPtr
upscaling_stagger
=
std
::
make_shared
<
CarboUspcaling
>
(
the_database
,
units_set
);
// Solver
// ------
// This section initialize the reactive transport solver
ReactiveTransportSolver
solver
(
transport_stagger
,
equilibrium_stagger
,
upscaling_stagger
);
set_reactive_solver_options
(
solver
.
get_options
());
transport_stagger
->
initialize
(
variables
);
equilibrium_stagger
->
initialize
(
variables
);
upscaling_stagger
->
initialize
(
variables
);
// Some basic information about the simulation
// First argument : the name of the simulation
// Second argument : the output function will be called every 3600s
SimulationInformation
simul_info
(
"carbo"
,
3600
);
simul_info
.
output_prefix
=
"out_carbo_"
;
// Runner
// -------
// Create the runner : loop through the timestep
ReactiveTransportRunner
runner
(
solver
,
0.1
,
100.0
,
simul_info
);
// Create the output class
OutputPolicy
output_policy
(
simul_info
,
the_mesh
,
the_database
,
specmicp_options
.
units_set
);
// and bind it to the runner
runner
.
set_output_policy
(
std
::
bind
(
&
OutputPolicy
::
output
,
&
output_policy
,
std
::
placeholders
::
_1
,
std
::
placeholders
::
_2
));
// Solve the problem
runner
.
run_until
(
1.0
*
24
*
3600
,
cast_var_to_base
(
variables
));
if
(
not
the_database
->
is_valid
())
throw
std
::
runtime_error
(
"The database has been changed during the computation"
);
// output information at the end of the timestep
io
::
print_minerals_profile
(
the_database
,
variables
,
the_mesh
,
simul_info
.
complete_filepath
(
"profiles_end"
,
"dat"
));
io
::
print_components_total_aqueous_concentration
(
the_database
,
variables
,
the_mesh
,
specmicp_options
.
units_set
,
simul_info
.
complete_filepath
(
"c_profiles_end"
,
"dat"
));
io
::
print_components_total_solid_concentration
(
the_database
,
variables
,
the_mesh
,
specmicp_options
.
units_set
,
simul_info
.
complete_filepath
(
"s_profiles_end"
,
"dat"
));
}
AdimensionalSystemSolverOptions
get_specmicp_options
()
{
AdimensionalSystemSolverOptions
options
;
options
.
solver_options
.
maxiter_maxstep
=
100
;
options
.
solver_options
.
maxstep
=
50.0
;
options
.
solver_options
.
threshold_cycling_linesearch
=
1e-6
;
options
.
solver_options
.
set_tolerance
(
1e-8
,
1e-14
);
options
.
solver_options
.
disable_condition_check
();
options
.
solver_options
.
disable_descent_direction
();
options
.
solver_options
.
enable_non_monotone_linesearch
();
options
.
solver_options
.
enable_scaling
();
options
.
system_options
.
non_ideality_tolerance
=
1e-14
;
options
.
system_options
.
cutoff_total_concentration
=
1e-10
;
options
.
units_set
.
length
=
units
::
LengthUnit
::
centimeter
;
return
options
;
}
void
set_transport_options
(
dfpmsolver
::
ParabolicDriverOptions
&
transport_options
)
{
transport_options
.
maximum_iterations
=
450
;
transport_options
.
quasi_newton
=
3
;
transport_options
.
residuals_tolerance
=
1e-8
;
transport_options
.
step_tolerance
=
1e-14
;
transport_options
.
absolute_tolerance
=
1e-18
;
transport_options
.
alpha
=
1.0
;
//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options
.
sparse_solver
=
SparseSolver
::
SparseLU
;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options
.
threshold_stationary_point
=
1e-4
;
}
void
set_reactive_solver_options
(
reactmicp
::
solver
::
ReactiveTransportOptions
&
reactive_options
)
{
reactive_options
.
maximum_iterations
=
50
;
//500;
reactive_options
.
residuals_tolerance
=
1e-2
;
reactive_options
.
good_enough_tolerance
=
1.0
;
reactive_options
.
absolute_residuals_tolerance
=
1e-14
;
reactive_options
.
implicit_upscaling
=
true
;
}
RawDatabasePtr
get_database
(
const
std
::
string
&
path_to_database
)
{
database
::
Database
the_database
(
path_to_database
);
std
::
map
<
std
::
string
,
std
::
string
>
swapping
({
{
"H[+]"
,
"HO[-]"
},
{
"Si(OH)4"
,
"SiO(OH)3[-]"
},
{
"Al[3+]"
,
"Al(OH)4[-]"
}
});
the_database
.
swap_components
(
swapping
);
the_database
.
remove_gas_phases
();
return
the_database
.
get_database
();
}
AdimensionalSystemConstraints
get_cement_initial_constraints
(
RawDatabasePtr
the_database
,
const
units
::
UnitsSet
&
units_set
)
{
// Vector init_min_publi(5);
// init_min_publi << 12.9, 1.62, 0.605, 0.131, 0.001;
// Vector init_min;
// init_min = initial_compo(init_min_publi);
// std::cout << "Init min \n ---- \n" << init_min << "\n ---- \n" << std::endl;
Vector
oxyde_compo
(
4
);
oxyde_compo
<<
66.98
,
21.66
,
7.19
,
2.96
;
Vector
species_compo
;
compo_from_oxyde
(
oxyde_compo
,
species_compo
);
scalar_t
mult
=
1e-6
;
//scalar_t mult = 6.73153e-4;
//scalar_t mult = 0.001296315;
//init_min *= 6.71146e-4;
//species_compo *= 0.947e-2; // poro 0.47
//species_compo *= 1.08e-2; // poro 0.4
//species_compo *= 1.25e-2; // poro 0.3
species_compo
*=
1.518e-2
;
Formulation
formulation
;
formulation
.
mass_solution
=
1.0
;
formulation
.
amount_minerals
=
{
{
"C3S"
,
species_compo
(
0
)},
{
"C2S"
,
species_compo
(
1
)},
{
"C3A"
,
species_compo
(
2
)},
{
"Gypsum"
,
species_compo
(
3
)},
{
"Calcite"
,
species_compo
(
0
)
*
1e-3
}
// {"Portlandite", init_min(0)},
// {"CSH,jennite", init_min(1)},
// {"Monosulfoaluminate", init_min(2)},
// {"Ettringite", init_min(3)},
// {"Calcite", init_min(4)}
};
formulation
.
extra_components_to_keep
=
{
"HCO3[-]"
};
formulation
.
concentration_aqueous
=
{
{
"Na[+]"
,
mult
*
1.0298
},
{
"K[+]"
,
mult
*
0.0801
},
{
"Cl[-]"
,
mult
*
0.0001
},
{
"HO[-]"
,
mult
*
1.8298
}
},
formulation
.
minerals_to_keep
=
{
"Portlandite"
,
"CSH,jennite"
,
"CSH,tobermorite"
,
"SiO2,am"
,
"Calcite"
,
"Al(OH)3,am"
,
"C3AH6"
,
"C4AH13"
,
"C2AH8"
,
"CAH10"
,
"Monosulfoaluminate"
,
"Tricarboaluminate"
,
"Monocarboaluminate"
,
"Hemicarboaluminate"
,
"Straetlingite"
,
"Gypsum"
,
"Ettringite"
,
"Thaumasite"
};
Dissolver
dissolver
(
the_database
);
dissolver
.
set_units
(
units_set
);
AdimensionalSystemConstraints
constraints
;
constraints
.
set_charge_keeper
(
1
);
constraints
.
total_concentrations
=
dissolver
.
dissolve
(
formulation
);
std
::
cout
<<
constraints
.
total_concentrations
<<
std
::
endl
;
constraints
.
set_saturated_system
();
return
constraints
;
}
AdimensionalSystemConstraints
get_bc_initial_constraints
(
RawDatabasePtr
the_database
,
const
units
::
UnitsSet
&
units_set
)
{
const
scalar_t
rho_w
=
laws
::
density_water
(
units
::
celsius
(
25.0
),
units_set
.
length
,
units_set
.
mass
);
const
scalar_t
nacl
=
0.3
*
rho_w
;
const
scalar_t
kcl
=
0.28
*
rho_w
;
Vector
total_concentrations
;
total_concentrations
.
setZero
(
the_database
->
nb_component
());
total_concentrations
(
the_database
->
get_id_component
(
"K[+]"
))
=
kcl
;
total_concentrations
(
the_database
->
get_id_component
(
"Cl[-]"
))
=
nacl
+
kcl
;
total_concentrations
(
the_database
->
get_id_component
(
"Na[+]"
))
=
nacl
;
total_concentrations
(
the_database
->
get_id_component
(
"SiO(OH)3[-]"
))
=
0.0001
*
rho_w
;
AdimensionalSystemConstraints
constraints
;
constraints
.
add_fixed_activity_component
(
the_database
->
get_id_component
(
"HO[-]"
),
-
10.3
);
constraints
.
set_charge_keeper
(
the_database
->
get_id_component
(
"HCO3[-]"
));
constraints
.
total_concentrations
=
total_concentrations
;
constraints
.
set_saturated_system
();
constraints
.
disable_surface_model
();
constraints
.
fixed_fugacity_cs
=
{};
return
constraints
;
}
AdimensionalSystemSolution
get_bc_initial_compo
(
RawDatabasePtr
the_database
,
const
AdimensionalSystemConstraints
&
constraints
,
const
AdimensionalSystemSolverOptions
&
options
)
{
Vector
variables
;
AdimensionalSystemSolver
solver
(
the_database
,
constraints
,
options
);
micpsolver
::
MiCPPerformance
perf
=
solver
.
solve
(
variables
,
true
);
if
(
perf
.
return_code
<=
micpsolver
::
MiCPSolverReturnCode
::
NotConvergedYet
)
throw
std
::
runtime_error
(
"Failed to solve initial composition"
);
std
::
cout
<<
variables
<<
std
::
endl
;
return
solver
.
get_raw_solution
(
variables
);
}
//AdimensionalSystemSolution get_cement_initial_compo(
// RawDatabasePtr the_database,
// const AdimensionalSystemConstraints& constraints,
// const AdimensionalSystemSolverOptions& options
// )
//{
// AdimensionalSystemSolver solver(the_database,
// constraints,
// options);
// Vector variables;
// solver.initialise_variables(
// variables,
// 0.4,
// {
// {"HO[-]", -1.2},
// {"Al(OH)4[-]", -4},
// {"HCO3[-]", -5},
// {"Ca[2+]", -1.6}
// });
// micpsolver::MiCPPerformance perf = solver.solve(variables);
// if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
// throw std::runtime_error("Failed to solve initial composition");
// std::cout << variables << std::endl;
// return solver.get_raw_solution(variables);
//}
void
compo_from_oxyde
(
Vector
&
compo_oxyde
,
Vector
&
compo_species
)
{
constexpr
double
M_CaO
=
56.08
;
constexpr
double
M_SiO2
=
60.09
;
constexpr
double
M_Al2O3
=
101.96
;
constexpr
double
M_SO3
=
80.06
;
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
::
ColPivHouseholderQR
<
Eigen
::
MatrixXd
>
solver
(
compo_in_oxyde
);
compo_species
=
solver
.
solve
(
compo_oxyde
);
}
Vector
initial_compo
(
const
Vector
&
publi
)
{
Matrix
publi_stochio
(
5
,
5
);
publi_stochio
<<
1
,
0
,
0
,
0
,
0
,
9
,
6
,
0
,
0
,
0
,
3
,
0
,
1
,
1
,
0
,
3
,
0
,
1
,
3
,
0
,
0
,
0
,
0
,
0
,
1
;
Matrix
cemdata_stochio
(
5
,
5
);
cemdata_stochio
<<
1
,
0
,
0
,
0
,
0
,
1.0
/
0.6
,
1.0
,
0
,
0
,
0
,
3
,
0
,
1
,
1
,
0
,
3
,
0
,
1
,
3
,
0
,
0
,
0
,
0
,
0
,
1
;
Vector
oxyde
;
Eigen
::
ColPivHouseholderQR
<
Matrix
>
solver
(
publi_stochio
);
oxyde
=
solver
.
solve
(
publi
);
Vector
for_cemdata
=
cemdata_stochio
*
oxyde
;
return
for_cemdata
;
}
mesh
::
Mesh1DPtr
get_mesh
(
scalar_t
dx
)
{
const
scalar_t
total_length
=
0.75
/
2.0
+
dx
;
//cm
const
scalar_t
height
=
20.0
;
index_t
nb_nodes
=
total_length
/
dx
+
1
;
return
mesh
::
uniform_axisymmetric_mesh1d
(
nb_nodes
,
total_length
,
dx
,
height
);
}
mesh
::
Mesh1DPtr
get_mesh
()
{
Vector
radius
(
71
);
const
scalar_t
height
=
20
;
// radius <<
// 3.751,
// 3.750,
// 3.749,
// 3.748,
// 3.747,
// 3.745,
// 3.740,
// 3.735,
// 3.730,
// 3.720,
// 3.710,
// 3.700,
// 3.675,
// 3.650,
// 3.625,
// 3.600,
// 3.575,
// 3.550,
// 3.525,
// 3.500,
// 3.475,
// 3.450,
// 3.425,
// 3.400,
// 3.375,
// 3.350,
// 3.325,
// 3.300,
// 3.250,
// 3.200,
// 3.150,
// 3.100,
// 3.050,
// 3.000,
// 2.975,
// 2.950,
// 2.925,
// 2.900,
// 2.875,
// 2.850,
// 2.825,
// 2.800,
// 2.775,
// 2.750,
// 2.700,
// 2.650;
radius
<<
3.751
,
3.750
,
3.745
,
3.740
,
3.735
,
3.730
,
3.720
,
3.710
,
3.700
,
3.690
,
3.680
,
3.670
,
3.660
,
3.645
,
3.630
,
3.615
,
3.600
,
3.580
,
3.560
,
3.540
,
3.520
,
3.500
,
3.475
,
3.450
,
3.425
,
3.400
,
3.375
,
3.350
,
3.325
,
3.300
,
3.275
,
3.250
,
3.225
,
3.200
,
3.175
,
3.150
,
3.125
,
3.100
,
3.050
,
3.025
,
3.000
,
2.950
,
2.900
,
2.850
,
2.800
,
2.750
,
2.700
,
2.650
,
2.600
,
2.550
,
2.500
,
2.450
,
2.400
,
2.350
,
2.300
,
2.250
,
2.200
,
2.150
,
2.100
,
2.050
,
2.000
,
1.950
,
1.900
,
1.850
,
1.800
,
1.750
,
1.700
,
1.650
,
1.600
,
1.550
,
1.500
;
radius
=
radius
/
10
;
return
mesh
::
axisymmetric_mesh1d
(
radius
,
height
);
}
Event Timeline
Log In to Comment