Page MenuHomec4science
Files F84187146

adimensional_system_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Sep 21, 06:05

adimensional_system_solver.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp/adimensional/adimensional_system.hpp"
#include "specmicp_common/micpsolver/micpsolver.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/io/adimensional_system_solution_saver.hpp"
#include "specmicp/io/adimensional_system_solution_reader.hpp"
#include "specmicp_database/database.hpp"
#include <iostream>
specmicp::RawDatabasePtr get_test_simple_database()
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
});
thedatabase.swap_components(swapping);
std::vector<std::string> to_keep = {"HO[-]", "Ca[2+]"};
thedatabase.keep_only_components(to_keep);
thedatabase.remove_half_cell_reactions(std::vector<std::string>({"H2O", "HO[-]",})) ;
return thedatabase.get_database();
}
using namespace specmicp;
TEST_CASE("Solving adimensional system", "[specmicp, MiCP, program, adimensional, solver]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Error;
specmicp::RawDatabasePtr thedatabase = get_test_simple_database();
auto id_h2o = database::DataContainer::water_index();
auto id_oh = thedatabase->get_id_component("HO[-]");
auto id_ca = thedatabase->get_id_component("Ca[2+]");
auto id_ch = thedatabase->get_id_mineral("Portlandite");
SECTION("Solving simple case") {
Vector total_concentration = Vector::Zero(thedatabase->nb_component());
total_concentration(id_h2o) = 3.0e4;
total_concentration(id_oh) = 2.0e4;
total_concentration(id_ca) = 1.0e4;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 100;
options.maxiter_maxstep = 100;
options.disable_crashing();
options.enable_scaling();
options.disable_descent_direction();
std::shared_ptr<specmicp::AdimensionalSystem> system=
std::make_shared<specmicp::AdimensionalSystem>(thedatabase, constraints);
specmicp::micpsolver::MiCPSolver<specmicp::AdimensionalSystem> solver(system);
solver.set_options(options);
system->number_eq();
Vector x = Vector::Zero(system->total_variables());
x(system->ideq_w()) = 0.8;
x(system->ideq_paq(id_oh)) = -2.0;
x(system->ideq_paq(id_ca)) = -2.3;
system->compute_log_gamma(x);
system->set_secondary_concentration(x);
solver.solve(x);
CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.462142).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.820933).epsilon(1e-4));
CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4));
}
SECTION("Solving simple case - volume in cm^3") {
Vector total_concentration = Vector::Zero(thedatabase->nb_component());
total_concentration(id_h2o) = 0.03;
total_concentration(id_oh) = 0.02;
total_concentration(id_ca) = 0.01;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 10;
options.maxiter_maxstep = 100;
options.disable_crashing();
options.enable_scaling();
options.disable_descent_direction();
std::shared_ptr<specmicp::AdimensionalSystem> system=
std::make_shared<specmicp::AdimensionalSystem>(thedatabase, constraints);
units::UnitsSet unit_set;
unit_set.length = units::LengthUnit::centimeter;
system->set_units(unit_set);
system->number_eq();
specmicp::micpsolver::MiCPSolver<specmicp::AdimensionalSystem> solver(system);
solver.set_options(options);
Vector x = Vector::Zero(system->total_variables());
x(system->ideq_w()) = 0.8;
x(system->ideq_paq(id_oh)) = -2.0;
x(system->ideq_paq(id_ca)) = -2.3;
system->compute_log_gamma(x);
system->set_secondary_concentration(x);
solver.solve(x);
CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4));
CHECK(x(system->ideq_min(id_ch)) == Approx( 0.32966).epsilon(1e-4));
}
SECTION("Solving simple case - volume in cm^3, mmol") {
Vector total_concentration = Vector::Zero(thedatabase->nb_component());
total_concentration(id_h2o) = 1e3*0.03;
total_concentration(id_oh) = 1e3*0.02;
total_concentration(id_ca) = 1e3*0.01;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 10;
options.maxiter_maxstep = 100;
options.disable_crashing();
options.enable_scaling();
options.disable_descent_direction();
std::shared_ptr<specmicp::AdimensionalSystem> system=
std::make_shared<specmicp::AdimensionalSystem>(thedatabase, constraints);
units::UnitsSet unit_set;
unit_set.length = units::LengthUnit::centimeter;
unit_set.quantity = units::QuantityUnit::millimoles;
system->set_units(unit_set);
system->number_eq();
specmicp::micpsolver::MiCPSolver<specmicp::AdimensionalSystem> solver(system);
solver.set_options(options);
Vector x = Vector::Zero(system->total_variables());
x(system->ideq_w()) = 0.8;
x(system->ideq_paq(id_oh)) = -2.0;
x(system->ideq_paq(id_ca)) = -2.3;
system->compute_log_gamma(x);
system->set_secondary_concentration(x);
solver.solve(x);
CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4));
CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4));
CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4));
}
SECTION("Automatic solver") {
Vector total_concentration = Vector::Zero(thedatabase->nb_component());
total_concentration(id_h2o) = 0.03;
total_concentration(id_oh) = 0.02;
total_concentration(id_ca) = 0.01;
specmicp::Vector x;
specmicp::AdimensionalSystemConstraints constraints(total_concentration);
specmicp::AdimensionalSystemSolver solver(thedatabase, constraints);
solver.initialize_variables(x, 0.8, -2.0);
//x(solver.dof_surface()) = -HUGE_VAL;
solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter;
solver.get_options().solver_options.maxstep = 10.0;
solver.get_options().solver_options.maxiter_maxstep = 100;
solver.get_options().solver_options.use_crashing = false;
solver.get_options().solver_options.use_scaling = true;
solver.get_options().solver_options.disable_descent_direction();
solver.get_options().solver_options.factor_gradient_search_direction = 100;
auto perf = solver.solve(x);
REQUIRE(perf.return_code >= micpsolver::MiCPSolverReturnCode::Success);
auto unsafe_solution = solver.unsafe_get_raw_solution(x);
auto solution = solver.get_raw_solution(x);
AdimensionalSystemSolutionExtractor unsafe_extr(
unsafe_solution, thedatabase, solver.get_options().units_set);
AdimensionalSystemSolutionExtractor extr(
solution, thedatabase, solver.get_options().units_set);
// check that we have the good solution
CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4));
CHECK(extr.log_molality_component(id_oh) == Approx(-1.46214).epsilon(1e-4));
CHECK(extr.log_molality_component(id_ca) == Approx(-1.82093).epsilon(1e-4));
//CHECK(extr.free_surface_concentration() == -HUGE_VAL);
CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4));
// check that unsafe solution is ok
CHECK(extr.volume_fraction_water()
== Approx(unsafe_extr.volume_fraction_water()).epsilon(1e-8));
CHECK(extr.log_molality_component(id_oh)
== Approx(unsafe_extr.log_molality_component(id_oh)).epsilon(1e-8));
CHECK(extr.log_molality_component(id_ca)
== Approx(unsafe_extr.log_molality_component(id_ca)).epsilon(1e-8));
CHECK(extr.volume_fraction_mineral(id_ch)
== Approx(unsafe_extr.volume_fraction_mineral(id_ch)).epsilon(1e-8));
CHECK(solution.log_gamma.norm()
== Approx(unsafe_solution.log_gamma.norm()).epsilon(1e-8));
CHECK(solution.secondary_molalities.norm()
== Approx(unsafe_solution.secondary_molalities.norm()).epsilon(1e-8));
database::Database the_db(thedatabase);
the_db.save("test_db.yaml");
io::save_solution_yaml("test_solution.yaml", thedatabase, solution, "test_db.yaml");
RawDatabasePtr new_db;
auto solution2 = io::parse_solution_yaml("test_solution.yaml", new_db);
REQUIRE(new_db != nullptr);
REQUIRE(solution.main_variables(0) == Approx(solution2.main_variables(0)));
REQUIRE(solution.main_variables(id_oh) == Approx(solution2.main_variables(id_oh)));
REQUIRE(solution.main_variables(id_ca) == Approx(solution2.main_variables(id_ca)));
REQUIRE(extr.volume_fraction_mineral(id_ch) == Approx(solution2.main_variables(extr.dof_mineral(id_ch))));
REQUIRE(solution.log_gamma.norm() == Approx(solution2.log_gamma.norm()));
CHECK(new_db->get_hash() == thedatabase->get_hash());
}
}

Event Timeline

c4science · Help