Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/specmicp/adim/adimensional_system_problem_solver.cpp
d94792a4f65fmaster

Last Change
Raw File

adimensional_system_problem_solver.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/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 "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_database/database.hpp"
#include <iostream>
TEST_CASE("Solving adimensinal system with problem setting", "[specmicp, MiCP, program, adimensional, solver]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
SECTION("Solving problem with dissolver - simpler problem") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
formulation.mass_solution = 0.8;
formulation.amount_minerals = {{"Portlandite", 10.0}};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
total_concentrations *= 1e3;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = raw_data->get_id_component("HO[-]");
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().solver_options.maxstep = 50.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.factor_descent_condition = 1e-10;
solver.get_options().solver_options.factor_gradient_search_direction = 200;
specmicp::Vector x;
solver.initialise_variables(x, 0.8, -2);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set);
CHECK(extr.volume_fraction_water() == Approx(0.802367).epsilon(1e-4));
CHECK(extr.volume_fraction_mineral(raw_data->get_id_mineral("Portlandite")) == Approx(0.32947).epsilon(1e-2));
}
SECTION("Solving problem with dissolver") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
formulation.mass_solution = 0.8;
formulation.amount_minerals = {{"C3S", 0.7}, {"C2S", 0.3}};
formulation.extra_components_to_keep = {"HCO3[-]", };
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
total_concentrations *= 1e3;
total_concentrations(2) = 500;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = raw_data->get_id_component("HO[-]");
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().solver_options.maxstep = 20.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.factor_descent_condition = 1e-10;
solver.get_options().solver_options.factor_gradient_search_direction = 200;
specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral());
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
}
SECTION("Solving problem with reactant box") {
std::cout << "\n Solving with reactant box \n ------------------------ \n";
std::cout << "== with SI units == \n";
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
reactant_box.set_charge_keeper("HO[-]");
specmicp::AdimensionalSystemConstraints constraints = reactant_box.get_constraints(true);
specmicp::Vector total_concentrations = constraints.total_concentrations;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialise_variables(x, 0.5, -6);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(solution,raw_data, specmicp::units::SI_units);
std::cout << "== with non SI units : mmol/dm == \n";
specmicp::units::UnitsSet dmmol;
dmmol.length = specmicp::units::LengthUnit::decimeter;
dmmol.quantity = specmicp::units::QuantityUnit::millimoles;
reactant_box.set_units(dmmol);
specmicp::AdimensionalSystemConstraints constraints2 = reactant_box.get_constraints(false);
specmicp::Vector total_concentrations2 = constraints2.total_concentrations;
specmicp::AdimensionalSystemSolver solver2(raw_data, constraints2);
solver2.get_options().units_set = dmmol;
solver2.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::Vector x2;
solver2.initialise_variables(x2, 0.5, -6);
perf = solver2.solve(x2);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
specmicp::AdimensionalSystemSolution solution2 = solver2.get_raw_solution(x);
CHECK(solution.main_variables(0) == Approx(solution2.main_variables(0)).epsilon(1e-6));
CHECK(solution.main_variables(2) == Approx(solution2.main_variables(2)).epsilon(1e-6));
CHECK(solution.main_variables(5) == Approx(solution2.main_variables(5)).epsilon(1e-6));
CHECK(solution.ionic_strength == Approx(solution2.ionic_strength).epsilon(1e-6));
specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, dmmol);
for (auto idc: raw_data->range_component()) {
if (idc == specmicp::database::electron_id) continue;
// std::cout << "Id " << idc << " - " << raw_data->get_label_component(idc) << std::endl;
// std::cout << total_concentrations(idc) << " - "
// << total_concentrations2(idc) << " - "
// << extr.total_concentration(idc) << " ~ "
// << extr.volume_fraction_gas_phase() * extr.total_gaseous_concentration(idc) << " # "
// << extr.fugacity_gas(0) << " - "
// << extr2.total_concentration(idc)
// << std::endl;
CHECK(extr.total_concentration(idc) == Approx(total_concentrations(idc)).epsilon(1e-8));
CHECK(extr2.total_concentration(idc) == Approx(total_concentrations2(idc)).epsilon(1e-8));
CHECK(extr.total_concentration(idc) == Approx(extr2.total_concentration(idc)).epsilon(1e-8));
}
std::cout << "== with non SI units : mmol/cm == \n";
specmicp::units::UnitsSet cmmol;
cmmol.length = specmicp::units::LengthUnit::centimeter;
cmmol.quantity = specmicp::units::QuantityUnit::millimoles;
reactant_box.set_units(cmmol);
specmicp::AdimensionalSystemConstraints constraints3 = reactant_box.get_constraints(false);
specmicp::Vector total_concentrations3 = constraints3.total_concentrations;
specmicp::AdimensionalSystemSolver solver3(raw_data, constraints3);
solver3.get_options().units_set = cmmol;
solver3.get_options().solver_options.set_tolerance(1e-8, 1e-12);
specmicp::Vector x3;
solver3.initialise_variables(x3, 0.5, -6);
perf = solver3.solve(x3);
for (auto idc: raw_data->range_component()) {
if (idc == specmicp::database::electron_id) continue;
CHECK(total_concentrations(idc) == Approx(1e3*total_concentrations3(idc)).epsilon(1e-8));
}
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << "perf nb it : " << perf.nb_iterations << std::endl;
specmicp::AdimensionalSystemSolution solution3 = solver3.get_raw_solution(x);
CHECK(solution.main_variables(0) == Approx(solution3.main_variables(0)).epsilon(1e-6));
CHECK(solution.main_variables(2) == Approx(solution3.main_variables(2)).epsilon(1e-6));
CHECK(solution.main_variables(5) == Approx(solution3.main_variables(5)).epsilon(1e-6));
CHECK(solution.ionic_strength == Approx(solution3.ionic_strength).epsilon(1e-6));
specmicp::AdimensionalSystemSolutionExtractor extr3(solution3, raw_data, cmmol);
for (auto idc: raw_data->range_component()) {
CHECK(extr3.total_concentration(idc) == Approx(total_concentrations3(idc)).epsilon(1e-8));
}
}
SECTION("Reactant box fixed activity") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 15, "mol/kg");
reactant_box.set_charge_keeper("Cl[-]");
reactant_box.add_fixed_activity_component("H[+]", 1e-8);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialise_variables(x, 0.8, -6);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.pH() == Approx(8));
}
SECTION("Reactant box fixed molality") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_aqueous_species("CO2", 15, "mol/kg");
reactant_box.set_charge_keeper("Cl[-]");
reactant_box.add_fixed_molality_component("H[+]", 1e-8);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialise_variables(x, 0.8, -6);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.molality_component(raw_data->get_id_component("H[+]")) == Approx(1e-8));
}
SECTION("Reactant box fixed fugacity") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"Si(OH)4", "SiO(OH)3[-]"}
});
thedatabase.swap_components(swapping);
//thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_solid_phase("C3S", 350, "L");
reactant_box.add_solid_phase("C2S", 50, "L");
reactant_box.add_fixed_fugacity_gas("CO2(g)", "HCO3[-]", 1e-3);
reactant_box.set_charge_keeper("H[+]");
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialise_variables(x, 0.8, -4);
solver.get_options().solver_options.set_tolerance(1e-6, 1e-8);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units);
CHECK(extr.fugacity_gas(raw_data->get_id_gas("CO2(g)")) == Approx(1e-3));
}
}
c4science · Help