Homec4science
Diffusion SpecMiCP / ReactMiCP (kinetics) Browse

specmicp/tests/specmicp/adim/adimensional_system_carboalu.cpp
8c94dee2ed89kinetics

Last Change
Raw File

adimensional_system_carboalu.cpp
View Options

#include <iostream>
#include "catch.hpp"
#include "utils/log.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
TEST_CASE("Carboalu - using adimensional system ", "[Adimensional, Carbonation, carboalu]")
{
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
SECTION("Carboalu")
{
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);
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 5e3;
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.extra_components_to_keep = {"HCO3[-]", };
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::index_t id_h2o = 0;
specmicp::index_t id_ho = 1;
specmicp::index_t id_co2 = 3;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = 1;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 50;
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-10;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral);
x.setZero();
x(0) = 0.5;
x.segment(1, raw_data->nb_component-1).setConstant(-6.0);
specmicp::scalar_t dh2co3 = mult*0.01;
for (specmicp::index_t i: raw_data->range_component())
{
std::cout << raw_data->labels_basis[i] << std::endl;
}
specmicp::uindex_t tot_iter = 0;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
for (int k=0; k<250; ++k)
{
constraints.total_concentrations(id_h2o) += dh2co3;
constraints.total_concentrations(id_ho) -= dh2co3;
constraints.total_concentrations(id_co2) += dh2co3;
solver = specmicp::AdimensionalSystemSolver (raw_data, constraints, solution, options);
//std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl;
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
//std::cout << x << std::endl;
REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
//std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
solution = solver.get_raw_solution(x);
//specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
//std::cout << solution.main_variables(0) << std::endl;
//std::cout << 14+solution.main_variables(1) << std::endl;
}
std::cout << tot_iter << std::endl;
}
}
c4science · Help