Homec4science
Diffusion SpecMiCP / ReactMiCP (oxydoreduc) Browse

specmicp/examples/reactmicp/equilibrium_curve.cpp
ca6efa1a4157oxydoreduc

Last Change
Raw File

equilibrium_curve.cpp
View Options

#include <iostream>
#include "utils/log.hpp"
#include "reactmicp/equilibrium_curve/chemistry.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
void test_chemistry() {
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 7e3;
specmicp::scalar_t m_c3s = mult*0.7;
specmicp::scalar_t m_c2s = mult*0.3;
specmicp::scalar_t wc = 0.5;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
std::cout << spec_solver.get_equilibrium_curve(500, 1.0) << std::endl;
}
void test_chemistry_with_al()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 6e3;
specmicp::scalar_t m_c3s = mult*0.6;
specmicp::scalar_t m_c2s = mult*0.2;
specmicp::scalar_t m_c3a = mult*0.10;
specmicp::scalar_t m_gypsum = mult*0.10;
specmicp::scalar_t wc = 0.8;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.02)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
{"C3A", m_c3a},
{"Gypsum", m_gypsum}
};
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Al(OH)3,am",
"Monosulfoaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
};
for (specmicp::index_t component: raw_data->range_component())
{
std::cout << raw_data->labels_basis[component] << std::endl;
}
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 20.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
std::cout << spec_solver.get_equilibrium_curve(10.0, -750.0) << std::endl;
}
int main()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
//test_chemistry();
test_chemistry_with_al();
}
c4science · Help