diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 4d235cc..a1619b9 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,137 +1,142 @@
###### Installation #######################
set(EXAMPLES_INSTALL_DIR ${SHARE_INSTALL_DIR}/examples)
# SpecMiCP
# ========
set (SPECMICP_EXAMPLES
specmicp/adimensional/thermocarbo.cpp
specmicp/adimensional/equilibrium_curve.cpp
specmicp/adimensional/momas_thermo.cpp
specmicp/adimensional/carbofe.cpp
)
install(FILES ${SPECMICP_EXAMPLES}
DESTINATION ${EXAMPLES_INSTALL_DIR}/specmicp
)
# ReactMiCP
# =========
set (REACTMICP_SATURATED_EXAMPLES
reactmicp/equilibrium_curve.cpp
reactmicp/saturated_react/react_leaching.cpp
reactmicp/saturated_react/carbonation.cpp
+ reactmicp/saturated_react/carbonation.yaml # the configuration file
reactmicp/saturated_react/carbonationfe.cpp
reactmicp/saturated_react/leaching_kinetic.cpp
)
install(FILES ${REACTMICP_SATURATED_EXAMPLES}
DESTINATION ${EXAMPLES_INSTALL_DIR}/reactmicp/saturated
)
+file(COPY reactmicp/saturated_react/carbonation.yaml
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
+)
+
###### Build (optional) ####################
if (NOT SPECMICP_BUILD_EXAMPLE)
set(EXAMPLE_IS_OPTIONAL EXCLUDE_FROM_ALL)
else()
set(EXAMPLE_IS_OPTIONAL )
endif()
# SpecMiCP
# ========
# thermocarbo
add_executable(ex_adim_thermocarbo
${EXAMPLE_IS_OPTIONAL}
specmicp/adimensional/thermocarbo.cpp
)
target_link_libraries(ex_adim_thermocarbo ${SPECMICP_LIBS})
# carbofe
add_executable(ex_adim_carbofe
${EXAMPLE_IS_OPTIONAL}
specmicp/adimensional/carbofe.cpp
)
target_link_libraries(ex_adim_carbofe ${SPECMICP_LIBS})
# Momas thermodynamic example
add_executable(ex_adim_momas
${EXAMPLE_IS_OPTIONAL}
specmicp/adimensional/momas_thermo.cpp
)
target_link_libraries(ex_adim_momas ${SPECMICP_LIBS})
# equilibrium curve
add_executable(ex_adim_equilibriumcurve
${EXAMPLE_IS_OPTIONAL}
specmicp/adimensional/equilibrium_curve.cpp
)
target_link_libraries(ex_adim_equilibriumcurve ${SPECMICP_LIBS})
# ReactMiCP
# =========
# EquilibriumCurve
add_executable(ex_reactmicp_equilibriumcurve
${EXAMPLE_IS_OPTIONAL}
reactmicp/equilibrium_curve.cpp
)
target_link_libraries(ex_reactmicp_equilibriumcurve ${REACTMICP_LIBS})
# Leaching
add_executable(ex_reactmicp_leaching
${EXAMPLE_IS_OPTIONAL}
reactmicp/saturated_react/react_leaching.cpp
)
target_link_libraries(ex_reactmicp_leaching ${REACTMICP_LIBS})
# Momas benchmark
add_executable(ex_reactmicp_momas
${EXAMPLE_IS_OPTIONAL}
reactmicp/saturated_react/momas_benchmark.cpp
)
target_link_libraries(ex_reactmicp_momas ${REACTMICP_LIBS})
# Carbonation
add_executable(ex_reactmicp_carbo
${EXAMPLE_IS_OPTIONAL}
reactmicp/saturated_react/carbonation.cpp
)
target_link_libraries(ex_reactmicp_carbo ${REACTMICP_LIBS})
add_executable(ex_reactmicp_carbo_static
EXCLUDE_FROM_ALL
reactmicp/saturated_react/carbonation.cpp
)
set_source_files_properties(reactmicp/saturated_react/carbonation.cpp PROPERTIES COMPILE_FLAGS -flto)
target_link_libraries(ex_reactmicp_carbo_static ${REACTMICP_STATIC_LIBS})
# Carbonation with Fe
add_executable(ex_reactmicp_carbofe
${EXAMPLE_IS_OPTIONAL}
reactmicp/saturated_react/carbonationfe.cpp
)
target_link_libraries(ex_reactmicp_carbofe ${REACTMICP_LIBS})
# Leaching with hydration
add_executable(ex_reactmicp_leaching_kinetic
${EXAMPLE_IS_OPTIONAL}
reactmicp/saturated_react/leaching_kinetic.cpp
)
target_link_libraries(ex_reactmicp_leaching_kinetic ${REACTMICP_STATIC_LIBS})
diff --git a/examples/reactmicp/saturated_react/carbonation.cpp b/examples/reactmicp/saturated_react/carbonation.cpp
index a39656a..09333e1 100644
--- a/examples/reactmicp/saturated_react/carbonation.cpp
+++ b/examples/reactmicp/saturated_react/carbonation.cpp
@@ -1,635 +1,510 @@
// ======================================================
// Leaching of a cement paste in CO2-saturated water
// =======================================================
// This file can be used to learn how to use ReactMiCP
// Include ReactMiCP
#include "reactmicp.hpp"
-#include "specmicp/adimensional/adimensional_system_solution_saver.hpp"
-#include "reactmicp/systems/saturated_react/diffusive_upscaling_stagger.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);
+// a uniform mesh (cf config)
// 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);
-
// ===================
// 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;
+
+ // Configuration
+
+ auto config = specmicp::io::parse_config_file("carbonation.yaml");
+
// The database
// -------------
// Obtain the database
- RawDatabasePtr the_database = get_database("../data/cemdata.js");
+ RawDatabasePtr the_database = specmicp::io::configure_database(config);
// SpecMiCP options
// ----------------
// Set the options
- AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
+ AdimensionalSystemSolverOptions specmicp_options;
+ specmicp::io::configure_specmicp_options(specmicp_options, config);
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
+ mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh(config);
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
// The equilibrium stagger
// - - - - - - - - - - - -
std::shared_ptr<EquilibriumStagger> equilibrium_stagger =
std::make_shared<satdiff::EquilibriumStagger>(list_constraints, index_constraints, specmicp_options);
// The transport stagger
// - - - - - - - - - - -
std::vector<index_t> list_fixed_component {0, 1, the_database->get_id_component("Si(OH)4")};
std::map<index_t, index_t> slave_nodes {};
std::shared_ptr<SaturatedTransportStagger> transport_stagger =
std::make_shared<satdiff::SaturatedTransportStagger>(
variables, list_fixed_nodes, slave_nodes, list_fixed_component
);
- set_transport_options(transport_stagger->options_solver());
+
+ specmicp::io::configure_transport_options(transport_stagger->options_solver(), config);
// The upscaling stagger
// - - - - - - - - - - -
CappedPowerLawParameters upscaling_param;
upscaling_param.d_eff_0 = 2.726e-12*1e4;
upscaling_param.cap = 1e10*1e6;
upscaling_param.exponent = 3.32;
upscaling_param.porosity_0 = 0.25;
upscaling_param.porosity_res = 0.02;
CappedPowerLaw upscaling_law(upscaling_param);
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<DiffusiveUpscalingStagger>(upscaling_law.get_law(),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());
+ specmicp::io::configure_reactmicp_options(solver.get_options(), config);
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_";
+ SimulationInformation simul_info = specmicp::io::configure_simulation_information(config);
// Runner
// -------
// Create the runner : loop through the timestep
- ReactiveTransportRunner runner(solver, 0.1, 100.0, simul_info);
+
+ specmicp::io::check_mandatory_yaml_node(config, "run", "__main__");
+ scalar_t min_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "minimum_dt", "run");
+ scalar_t max_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "maximum_dt", "run");
+ scalar_t total = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "total_execution_time", "run");
+ ReactiveTransportRunner runner(solver, min_dt, max_dt, simul_info);
// Output
// ------
io::OutputNodalVariables output_policy(the_database, the_mesh, specmicp_options.units_set);
+ io::configure_reactmicp_output(output_policy, simul_info, config);
- std::string suffix { "dat" };
-
- index_t id_hco3 = the_database->get_id_component("HCO3[-]");
- index_t id_ca = the_database->get_id_component("Ca[2+]");
- index_t id_calcite = the_database->get_id_mineral("Calcite");
- index_t id_si = the_database->get_id_component("Si(OH)4");
-
- output_policy.register_porosity(simul_info.complete_filepath("poro", suffix));
- output_policy.register_total_aqueous_concentration(id_hco3, simul_info.complete_filepath("c_hco3", suffix));
- output_policy.register_total_solid_concentration(id_ca, simul_info.complete_filepath("s_ca", suffix));
- output_policy.register_diffusion_coefficient(simul_info.complete_filepath("diffusivity", suffix));
- output_policy.register_pH(simul_info.complete_filepath("ph", suffix));
- output_policy.register_volume_fraction_mineral(id_calcite, simul_info.complete_filepath("calcite", suffix));
- output_policy.register_total_concentration(id_si, simul_info.complete_filepath("t_si", suffix));
// and bind it to the runner
runner.set_output_policy(output_policy.get_output_for_reactmicp());
// Solve the problem
// -----------------
- runner.run_until(5.0*24*3600, cast_var_to_base(variables));
+ runner.run_until(total, cast_var_to_base(variables));
if(not the_database->is_valid())
throw std::runtime_error("The database has been changed during the computation");
AdimensionalSystemSolutionSaver saver(the_database);
saver.save_solutions(simul_info.complete_filepath("solutions", "dat"), variables->equilibrium_solutions());
// output information at the end of the timestep
io::print_minerals_profile(the_database, variables, the_mesh, specmicp_options.units_set,
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"));
return EXIT_SUCCESS;
}
-
-AdimensionalSystemSolverOptions get_specmicp_options()
-{
- AdimensionalSystemSolverOptions options;
-
-
- options.solver_options.maxiter_maxstep = 100;
- options.solver_options.maxstep = 20.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 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;
//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.1e-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}
};
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(mic)",
- "C3AH6",
- "C3ASO_41H5_18",
- "C3ASO_84H4_32",
- "C4AH19",
- "C4AH13",
- "C2AH7_5",
- "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("Si(OH)4")) = 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);
}
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)
-{
- mesh::UniformAxisymmetricMesh1DGeometry geometry;
- geometry.dx = dx;
- geometry.radius = 0.75/2.0 + dx; //cm
- geometry.height = 20.0;
- geometry.nb_nodes = geometry.radius/dx + 1;
-
- return mesh::uniform_axisymmetric_mesh1d(geometry);
-}
-
mesh::Mesh1DPtr get_mesh()
{
Vector radius(71);
const scalar_t height = 20;
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);
}
diff --git a/examples/reactmicp/saturated_react/carbonation.yaml b/examples/reactmicp/saturated_react/carbonation.yaml
new file mode 100644
index 0000000..07707d8
--- /dev/null
+++ b/examples/reactmicp/saturated_react/carbonation.yaml
@@ -0,0 +1,111 @@
+# --------------------------------------------------#
+# Leaching of a cement paste in CO2-saturated water #
+# --------------------------------------------------#
+
+# This is an example of a configuration file for ReactMiCP
+
+# These options are used to buld the database
+database:
+ path: "../data/cemdata.js" # Path to the database
+ swap_components: # Set the bases
+ - {out: "H[+]", in: "HO[-]"}
+ - {out: "Al[3+]", in: "Al(OH)4[-]"}
+ remove_gas: true # remove the gas species
+ list_solids_to_keep: # the solids phases at equilibrium to keep
+ - "Portlandite"
+ - "CSH,jennite"
+ - "CSH,tobermorite"
+ - "SiO2(am)"
+ - "Calcite"
+ - "Al(OH)3(mic)"
+ - "C3AH6"
+ - "C3ASO_41H5_18"
+ - "C3ASO_84H4_32"
+ - "C4AH19"
+ - "C4AH13"
+ - "C2AH7_5"
+ - "CAH10"
+ - "Monosulfoaluminate"
+ - "Tricarboaluminate"
+ - "Monocarboaluminate"
+ - "Hemicarboaluminate"
+ # - "Straetlingite"
+ - "Gypsum"
+ - "Ettringite"
+ # - "Thaumasite"
+
+# The units
+units:
+ length: centimeter
+
+# These options are used to build the mesh
+mesh:
+ type: uniform_axisymmetric
+ uniform_axisymmetric:
+ dx: 0.005
+ radius: 0.378
+ height: 20.0
+
+# Options for SpecMiCP
+specmicp_options:
+ maximum_step_length: 20.0
+ maximum_step_max_iteration: 100
+ threshold_cycling_linesearch: 1e-6
+ residual_tolerance: 1e-8
+ step_tolerance: 1e-14
+ threshold_condition_check: -1
+ factor_descent_direction: -1
+ enable_nonmonotone_linesearch: true
+ enable_scaling: true
+ non_ideality_tolerance: 1e-14
+ cutoff_total_concentration: 1e-10
+
+# Options for the transport solver
+transport_options:
+ maximum_iteration: 450
+ quasi_newton: 3
+ residual_tolerance: 1e-8
+ absolute_tolerance: 1e-14
+ step_tolerance: 1e-14
+ sparse_solver: LU
+ # linesearch: strang
+ threshold_stationary: 1e-4
+
+# options for the reactive transport solver
+reactmicp_options:
+ residual_tolerance: 1e-2
+ absolute_tolerance: 1e-14
+ # step_tolerance:
+ good_enough_tolerance: 1.0
+ maximum_iteration: 50
+ implicit_upscaling: true
+
+simulation:
+ name: carbo
+ output_prefix: out_carbo_
+ print_iter_info: true
+ output_step: 3600
+
+run:
+ minimum_dt: 0.1
+ maximum_dt: 100.0
+ total_execution_time: 432000 #5.0*3600*24
+
+reactmicp_output:
+ porosity: true
+ diffusivity: true
+ total_aqueous_concentration:
+ - "HCO3[-]"
+ - "Ca[2+]"
+ total_solid_concentration:
+ - "Ca[2+]"
+ - "Al(OH)4[-]"
+ ph: true
+ volume_fraction_solid:
+ - "Calcite"
+ total_concentration:
+ - "Si(OH)4"
+
+
+
+