Homec4science
Diffusion SpecMiCP / ReactMiCP (chloride) Browse

specmicp/tests/specmicp/adim/adimensional_system_thermocarbo.cpp
88a299bc1607chloride

Last Change
Raw File

adimensional_system_thermocarbo.cpp
View Options

#include "catch.hpp"
#include <iostream>
#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"
#include "utils/timer.hpp"
#include "database/database.hpp"
TEST_CASE("thermocarbo - using adimensional system ", "[Adimensional, Thermocarbo]")
{
std::cerr.flush();
std::cout.flush();
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
SECTION("Thermocarbo")
{
std::cout << "\n-------------------\n Thermocarbo\n -------------------" << std::endl;
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 1.0;
formulation.mass_solution = mult*0.156;
formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}};
formulation.extra_components_to_keep = {"HCO3[-]", };
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = raw_data->get_id_component("HO[-]");
specmicp::index_t id_h2o = raw_data->water_index();
specmicp::index_t id_ho = raw_data->get_id_component("HO[-]");
specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]");
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 100.0;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.disable_descent_direction();
options.solver_options.set_tolerance(1e-8, 1e-12);
//options.solver_options.disable_non_monotone_linesearch();
options.solver_options.disable_condition_check();
options.solver_options.disable_crashing();
options.solver_options.disable_scaling();
specmicp::Timer tot_timer;
tot_timer.start();
specmicp::Vector x;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(x, 0.8, -4.0);
specmicp::scalar_t dh2co3 = 0.1;
specmicp::index_t nb_iter {0};
constexpr specmicp::index_t nb_step = 40;
specmicp::Vector solution(nb_step);
for (int k=0; k<nb_step; ++k)
{
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(static_cast<int>(perf.return_code) >=
static_cast<int>(specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet));
nb_iter += perf.nb_iterations;
constraints.total_concentrations(id_h2o) += mult*dh2co3;
constraints.total_concentrations(id_ho) -= mult*dh2co3;
constraints.total_concentrations(id_co2) += mult*dh2co3;
specmicp::AdimensionalSystemSolution adim_solution = solver.get_raw_solution(x);
solution(k) = 14+adim_solution.main_variables(id_ho);
}
tot_timer.stop();
std::cout << "Total time : " << tot_timer.elapsed_time() << "s (Ref 0.040s)" << std::endl;
std::cout << "Nb iter : " << nb_iter << " (Ref 1846)" << std::endl;
specmicp::Vector reference_solution(nb_step);
reference_solution << 12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.5194,
12.1374,
12.1374,
12.1374,
12.1374,
12.1374,
12.1374,
12.1374,
12.1374,
9.8575,
9.8575,
9.8575,
9.8575,
9.8575,
9.8575,
9.8575,
9.8575,
8.97554,
5.33701,
5.15035,
5.0435,
4.96886,
4.91174,
4.86562,
4.82705,
4.794,
4.76515,
4.73958,
4.71668,
4.69597;
for (int k=0;k<nb_step; ++k) {
CHECK(solution(k) == Approx(reference_solution(k)).epsilon(1e-3));
}
}
}
c4science · Help