Page MenuHomec4science
Files F84887656

adimensional_system_problem_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, Sep 25, 10:16

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/problem_solver/step_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 adimensional 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 immobilized immobile species") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]", "Al(OH)4[-]"},
{"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::units::UnitsSet dm_mol;
dm_mol.length = specmicp::units::LengthUnit::decimeter;
specmicp::ReactantBox reactant_box(raw_data, dm_mol);
reactant_box.set_solution(0.5, "L");
reactant_box.add_aqueous_species("Ca(OH)2", 0.1, "mol/kg");
reactant_box.add_aqueous_species("Na2SO4", 0.5, "mol/kg");
reactant_box.disable_immobile_species();
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);
solver.get_options().units_set = dm_mol;
solver.get_options().solver_options.set_tolerance(1e-8, 1e-10);
specmicp::Vector x;
solver.initialize_variables(x, 0.5, -6);
auto perf = solver.solve(x);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, dm_mol);
for (auto m: raw_data->range_mineral()) {
CHECK(extr.volume_fraction_mineral(m) == 0.0);
CHECK(extr.saturation_index(m) > 0.0);
}
}
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("Solving problem with reactant box and smart solver CEMDATA18") {
std::cout << "\n Solving with reactant box and smart solver \n ------------------------ \n";
specmicp::database::Database thedatabase(TEST_CEMDATA18_PATH);
std::map<std::string, std::string> swapping ({
{"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(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("Solving problem with reactant box and smart solver, fixed saturation") {
std::cout << "\n Solving with reactant box and smart solver, fixed saturation \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[-]");
reactant_box.set_fixed_saturation(0.8);
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);
auto sol = solver.get_solution();
auto extr = specmicp::AdimensionalSystemSolutionExtractor(sol, raw_data, reactant_box.get_units());
CHECK(extr.saturation_water() == Approx(0.8).epsilon(1e-6));
}
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));
}
SECTION("Reactant box - fixed saturation index") {
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units);
reactant_box.set_solution(1, "L/L");
reactant_box.add_fixed_saturation_index("Calcite", "HCO3[-]", 0.0);
reactant_box.add_fixed_activity_component("H[+]", 1e-4);
reactant_box.set_charge_keeper("Ca[2+]");
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.pH() == Approx(4.0));
CHECK(extr.saturation_index(raw_data->get_id_mineral("Calcite")) == Approx(0.0).margin(1e-6));
}
SECTION("Reactant box - mineral upper bound") {
std::cout << "\n Reactant box : box vi \n ------------------------ \n";
std::cout << "first solution" << std::endl;
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
std::map<std::string, std::string> swapping ({
//{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]", "Al(OH)4[-]"},
{"HCO3[-]", "CO2"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"});
thedatabase.minerals_keep_only({
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2(am)",
"Calcite",
"Al(OH)3(mic)",
"C3AH6",
"C3ASO_41H5_18",
"C3ASO_84H4_32",
"C4AH19",
"C4AH13",
"C2AH7_5",
"CAH10",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Gypsum",
"Ettringite"
});
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::units::UnitsSet units_set;
units_set.length = specmicp::units::LengthUnit::centimeter;
units_set.quantity = specmicp::units::QuantityUnit::millimoles;
specmicp::ReactantBox reactant_box(raw_data, units_set);
reactant_box.set_solution(0.25, "kg/L");
reactant_box.add_solid_phase("Portlandite", 4.46, "mol/L");
reactant_box.add_solid_phase("CSH,jennite", 4.59, "mol/L");
reactant_box.add_solid_phase("Monosulfoaluminate", 0.3, "mol/L");
reactant_box.add_solid_phase("Al(OH)3(mic)", 0.6, "mol/L");
reactant_box.add_solid_phase("Calcite", 0.2, "mol/L");
//reactant_box.add_aqueous_species("CO2", 0.5, "mol/L");
reactant_box.add_aqueous_species("KOH", 473.0, "mmol/L");
reactant_box.add_aqueous_species("NaOH", 52.0, "mmol/L");
reactant_box.set_fixed_saturation(0.7);
auto constraints = reactant_box.get_constraints(true);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints);
solver.get_options().units_set = units_set;
solver.get_options().system_options.non_ideality_tolerance = 1e-12;
specmicp::Vector x;
solver.initialize_variables(x, 0.7, -6);
// /solver.get_options().solver_options.set_tolerance(1e-6, 1e-8);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, units_set);
std::cout << "constrained solution" << std::endl;
specmicp::index_t id_ettr = raw_data->get_id_mineral("Ettringite");
specmicp::scalar_t max_vol_frac = extr.volume_fraction_mineral(id_ettr);
reactant_box.set_mineral_upper_bound("Ettringite", max_vol_frac);
reactant_box.add_solid_phase("Gypsum", 0.5, "mol/L");
constraints = reactant_box.get_constraints(false);
solver = specmicp::AdimensionalSystemSolver(raw_data, constraints);
solver.get_options().units_set = units_set;
perf = solver.solve(x, false);
REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
auto sol2 = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor extr2(sol2, raw_data, units_set);
CHECK(extr2.volume_fraction_mineral(id_ettr) <= max_vol_frac);
std::cout << extr2.volume_fraction_mineral(id_ettr) << " ?(<=) " << max_vol_frac << std::endl;
}
SECTION("Solving problem with stepper") {
std::cout << "\n Solving with stepper \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();
// initial constraints
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.set_charge_keeper("HO[-]");
reactant_box.add_component("CO2");
specmicp::uindex_t max_step = 20;
specmicp::StepProblem problem;
problem.set_initial_constraints(reactant_box, true);
specmicp::index_t id_co2 = raw_data->get_id_component("CO2");
specmicp::index_t id_ca = raw_data->get_id_component("Ca[2+]");
// The updates to apply at each step
specmicp::Vector update;
update.setZero(raw_data->nb_component());
update(id_co2) = (0.95*problem.get_constraints().total_concentrations(id_ca))/max_step;
problem.set_constraint_updater(specmicp::StepTotalConcentration(update));
problem.set_stop_condition_step(max_step);
specmicp::StepProblemRunner runner;
// output
std::vector<double> pHs;
pHs.reserve(max_step);
runner.set_step_output(
[&pHs](specmicp::uindex_t cnt,
const specmicp::AdimensionalSystemSolutionExtractor& extr) {
pHs.push_back(extr.pH());
return specmicp::StepProblemStatus::OK;
});
auto status = runner.run(problem);
REQUIRE(status == specmicp::StepProblemStatus::OK);
CHECK(pHs.size() == max_step);
for (auto it: pHs) {
std::cout << it << std::endl;
}
CHECK(pHs[pHs.size()-1] == Approx(9.84065).epsilon(1e-4));
}
}

Event Timeline

c4science · Help