Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120472332
adimensional_system_problem_solver.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Jul 4, 15:26
Size
19 KB
Mime Type
text/x-c
Expires
Sun, Jul 6, 15:26 (2 d)
Engine
blob
Format
Raw Data
Handle
27193778
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Log In to Comment