Homec4science
Diffusion SpecMiCP / ReactMiCP (kinetics) Browse

specmicp/tests/specmicp/adim/adimensional_system_solver.cpp
7dcd83852458kinetics

Last Change
Raw File

adimensional_system_solver.cpp
View Options

#include "catch.hpp"
#include <iostream>
#include "utils/log.hpp"
#include "specmicp/adimensional/adimensional_system.hpp"
#include "micpsolver/micpsolver.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
specmicp::RawDatabasePtr get_test_simple_database()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
});
thedatabase.swap_components(swapping);
std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]"};
thedatabase.keep_only_components(to_keep);
return thedatabase.get_database();
}
TEST_CASE("Solving adimensional system", "[specmicp, MiCP, program, adimensional, solver]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Error;
SECTION("Solving simple case") {
specmicp::RawDatabasePtr thedatabase = get_test_simple_database();
specmicp::Vector total_concentrations(thedatabase->nb_component);
total_concentrations << 3.0e4, 2e4 , 1e4;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 100;
options.maxiter_maxstep = 100;
options.use_crashing = false;
options.use_scaling = true;
options.penalization_factor = 0.8;
options.max_factorization_step = 4;
options.factor_descent_condition = 1e-10;
std::shared_ptr<specmicp::AdimensionalSystem> system=
std::make_shared<specmicp::AdimensionalSystem>(thedatabase, constraints);
specmicp::micpsolver::MiCPSolver<specmicp::AdimensionalSystem> solver(system);
solver.set_options(options);
specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral);
x << 0.8, -2.0, -2.3, 0.0;
system->compute_log_gamma(x);
system->set_secondary_concentration(x);
solver.solve(x);
REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4));
REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4));
REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4));
REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4));
}
SECTION("Solving simple case - volum in cm^3") {
specmicp::RawDatabasePtr thedatabase = get_test_simple_database();
specmicp::Vector total_concentrations(thedatabase->nb_component);
total_concentrations << 0.03, 0.02, 0.01;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 100;
options.maxiter_maxstep = 100;
options.use_crashing = false;
options.use_scaling = false;
options.penalization_factor = 0.8;
options.max_factorization_step = 4;
options.factor_descent_condition = 1e-10;
std::shared_ptr<specmicp::AdimensionalSystem> system=
std::make_shared<specmicp::AdimensionalSystem>(thedatabase, constraints);
system->length_unit() = specmicp::units::LengthUnit::centimeter;
specmicp::micpsolver::MiCPSolver<specmicp::AdimensionalSystem> solver(system);
solver.set_options(options);
specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral);
x << 0.8, -2.0, -2.3, 0.0;
system->compute_log_gamma(x);
system->set_secondary_concentration(x);
solver.solve(x);
REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4));
REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4));
REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4));
REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4));
}
SECTION("Automatic solver") {
specmicp::RawDatabasePtr thedatabase = get_test_simple_database();
specmicp::Vector total_concentrations(thedatabase->nb_component);
total_concentrations << 0.03, 0.02, 0.01;
specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral);
x << 0.8, -2.0, -2.3, 0.0;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
specmicp::AdimensionalSystemSolver solver(thedatabase, constraints);
solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter;
solver.get_options().solver_options.maxstep = 100.0;
solver.get_options().solver_options.maxiter_maxstep = 100;
solver.get_options().solver_options.use_crashing = false;
solver.get_options().solver_options.use_scaling = false;
solver.get_options().solver_options.factor_descent_condition = 1e-10;
solver.get_options().solver_options.factor_gradient_search_direction = 100;
solver.solve(x);
REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4));
REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4));
REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4));
REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4));
}
}
c4science · Help