Page MenuHomec4science
Files F62663263

adimensional_system_problem_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, May 14, 17:09

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/problem_solver/smart_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_database/database.hpp"
#include "specmicp_common/physics/laws.hpp"
#include <iostream>
TEST_CASE("Reactant Box", "[units],[init],[formulation]")
{
SECTION("Units setting - SI")
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
auto raw_data= thedatabase.get_database();
auto units_set = specmicp::units::SI_units;
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(500, "L");
reactant_box.add_aqueous_species("CaCl2", 5.0, "kg");
reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg");
auto vector = reactant_box.get_total_concentration(false);
auto mass_sol = 0.5 * specmicp::laws::density_water(units_set);
CHECK(vector(0) == Approx( mass_sol
/ raw_data->molar_mass_basis(0, units_set)
));
auto id_ca = raw_data->get_id_component("Ca[2+]");
auto id_na = raw_data->get_id_component("Na[+]");
auto id_cl = raw_data->get_id_component("Cl[-]");
auto id_cacl2 = raw_data->get_id_compound("CaCl2");
CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2)));
CHECK(vector(id_na) == Approx(0.5*mass_sol));
CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca)));
}
SECTION("Units setting - mol L")
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
auto raw_data= thedatabase.get_database();
auto units_set = specmicp::units::SI_units;
units_set.length = specmicp::units::LengthUnit::decimeter;
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(0.5, "L");
reactant_box.add_aqueous_species("CaCl2", 5.0, "kg");
reactant_box.add_aqueous_species("NaCl", 0.5, "mol/kg");
auto vector = reactant_box.get_total_concentration(false);
auto mass_sol = 0.5 * specmicp::laws::density_water(units_set);
CHECK(vector(0) == Approx( mass_sol
/ raw_data->molar_mass_basis(0, units_set)
));
auto id_ca = raw_data->get_id_component("Ca[2+]");
auto id_na = raw_data->get_id_component("Na[+]");
auto id_cl = raw_data->get_id_component("Cl[-]");
auto id_cacl2 = raw_data->get_id_compound("CaCl2");
CHECK(vector(id_ca) == Approx(5.0/raw_data->molar_mass_compound(id_cacl2)));
CHECK(vector(id_na) == Approx(0.5*mass_sol));
CHECK(vector(id_cl) == Approx(vector(id_na) + 2*vector(id_ca)));
}
}
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.initialize_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.initialize_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.initialize_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.initialize_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("Solving problem with reactant box and smart solver") {
std::cout << "\n Solving with reactant box and smart solver \n ------------------------ \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::SmartAdimSolver solver(raw_data, reactant_box);
solver.set_init_volfrac_water(0.5);
solver.set_init_molality(1e-6);
solver.solve();
REQUIRE(solver.is_solved() == true);
}
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(0.5, "L/L");
reactant_box.add_solid_phase("C3S", 0.350, "L/L");
reactant_box.add_solid_phase("C2S", 0.05, "L/L");
reactant_box.set_charge_keeper("HCO3[-]");
reactant_box.add_fixed_activity_component("H[+]", 1e-9);
auto constraints = reactant_box.get_constraints(true);
CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]"));
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.6, -4);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
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(9));
}
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[-]"},
// {"H[+]", "HO[-]"}
});
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(0.5, "L/L");
reactant_box.add_solid_phase("C3S", 0.35, "L/L");
reactant_box.add_solid_phase("C2S", 0.05, "L/L");
reactant_box.set_charge_keeper("HCO3[-]");
reactant_box.add_fixed_molality_component("H[+]", 1e-5);
auto constraints = reactant_box.get_constraints(true);
CHECK(constraints.charge_keeper == raw_data->get_id_component("HCO3[-]"));
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
specmicp::Vector x;
solver.initialize_variables(x, 0.6, -4);
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
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-5));
}
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.initialize_variables(x, 0.8, -6);
solver.get_options().solver_options.set_tolerance(1e-6, 1e-8);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
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));
}
}

Event Timeline

c4science · Help