diff --git a/tests/specmicp/adim/adimensional_system.cpp b/tests/specmicp/adim/adimensional_system.cpp index 18d0425..8dbbffc 100644 --- a/tests/specmicp/adim/adimensional_system.cpp +++ b/tests/specmicp/adim/adimensional_system.cpp @@ -1,147 +1,175 @@ #include "catch.hpp" #include "specmicp/adimensional/adimensional_system.hpp" +#include "database/database.hpp" #include specmicp::RawDatabasePtr get_test_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); - std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]"}; + std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); return thedatabase.get_database(); } +using namespace specmicp; TEST_CASE("Adimensional system", "[specmicp, MiCP, program, adimensional]") { - specmicp::RawDatabasePtr thedatabase = get_test_database(); + RawDatabasePtr thedatabase = get_test_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("initialization") { - specmicp::Vector total_concentration(thedatabase->nb_component); - total_concentration << 55.5, 2e-3, 1e-3; - specmicp::AdimensionalSystemConstraints constraints(total_concentration); - specmicp::AdimensionalSystem system(thedatabase, constraints); + Vector total_concentration = Vector::Zero(thedatabase->nb_component()); + total_concentration(id_h2o) = 55.5; + total_concentration(id_oh) = 2e-3; + total_concentration(id_ca) = 1e-3; + AdimensionalSystemConstraints constraints(total_concentration); + + AdimensionalSystem system(thedatabase, constraints); - REQUIRE(system.ideq_w() == 0); - REQUIRE(system.ideq_paq(1) == 1); - REQUIRE(system.ideq_paq(2) == 2); - REQUIRE(system.ideq_min(0) == 3); + REQUIRE(system.ideq_w() != no_equation); + REQUIRE(system.ideq_paq(id_oh) != no_equation); + REQUIRE(system.ideq_paq(id_ca) != no_equation); + REQUIRE(system.ideq_min(id_ch) != no_equation); - specmicp::Vector x(5); - x << 1.0, std::log10(2e-3), -3, 0.0, 0.0; + Vector x = Vector::Zero(system.total_dofs()); + x(id_h2o) = 1.0; + x(id_oh) = std::log10(2e-3); + x(id_ca) = -3; system.set_secondary_concentration(x); - specmicp::scalar_t molality_h = system.secondary_molality(0); + scalar_t molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); REQUIRE(molality_h == Approx(5e-12)); system.compute_log_gamma(x); molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); system.hook_start_iteration(x, 1.0); molality_h = system.secondary_molality(0); REQUIRE(std::isfinite(molality_h)); system.hook_start_iteration(x, 1.0); REQUIRE(molality_h == Approx(system.secondary_molality(0))); } SECTION("Residuals") { - specmicp::Vector total_concentration(thedatabase->nb_component); - total_concentration << 55.5, 2e-3, 1e-3; + Vector total_concentration = Vector::Zero(thedatabase->nb_component()); + total_concentration(id_h2o) = 55.5; + total_concentration(id_oh) = 2e-3; + total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystem system(thedatabase, constraints); - specmicp::Vector x(4); - x << 1.0, std::log10(2e-3), -3, 0.0; + Vector x = Vector::Zero(system.total_dofs()); + x(system.ideq_w()) = 1.0; + x(system.ideq_paq(id_oh)) = std::log10(2e-3); + x(system.ideq_paq(id_ca)) = -3; specmicp::scalar_t res = system.residual_water(x); REQUIRE(std::isfinite(res)); - REQUIRE(res == Approx((55.5 - system.density_water()/thedatabase->molar_mass_basis_si(0))/55.5)); + REQUIRE(res == Approx((55.5 - system.density_water()/thedatabase->molar_mass_basis_si(id_h2o))/55.5)); } SECTION("Equation numbering") { - specmicp::Vector total_concentration(thedatabase->nb_component); - total_concentration << 55.5, 2e-3, 1e-3; + Vector total_concentration = Vector::Zero(thedatabase->nb_component()); + total_concentration(id_h2o) = 55.5; + total_concentration(id_oh) = 2e-3; + total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.disable_conservation_water(); constraints.enable_surface_model(1.23456); constraints.total_concentrations(2) = 0; specmicp::AdimensionalSystem system(thedatabase, constraints); CHECK(system.water_equation_type() == specmicp::WaterEquationType::NoEquation); CHECK(system.surface_model() == specmicp::SurfaceEquationType::Equilibrium); CHECK(system.surface_total_concentration() == 1.23456); CHECK((int) system.aqueous_component_equation_type(2) == (int) specmicp::AqueousComponentEquationType::NoEquation); constraints = specmicp::AdimensionalSystemConstraints(total_concentration); constraints.enable_conservation_water(); constraints.disable_surface_model(); constraints.add_fixed_activity_component(2, -2.0); system = specmicp::AdimensionalSystem(thedatabase, constraints); CHECK(system.water_equation_type() == specmicp::WaterEquationType::MassConservation); CHECK(system.surface_model() == specmicp::SurfaceEquationType::NoEquation); CHECK(system.aqueous_component_equation_type(2) == specmicp::AqueousComponentEquationType::FixedActivity); CHECK(system.fixed_activity_bc(2) == -2.0); } SECTION("Jacobian") { // add a sorbed species, just to see if it works - thedatabase->nb_sorbed = 1; - thedatabase->nu_sorbed = specmicp::Matrix(1, thedatabase->nb_component+1); - thedatabase->nu_sorbed << 0, 2, 1, 1; - thedatabase->logk_sorbed = specmicp::Vector(1); - thedatabase->logk_sorbed << -1; - - - specmicp::Vector total_concentration(thedatabase->nb_component); - total_concentration << 55.5, 2e-3, 1e-3; + thedatabase->sorbed.resize(1); + specmicp::database::SorbedValues values; + values.label = "SorbedSpecies"; + values.logk = -1; + values.sorption_site_occupied = 1; + thedatabase->sorbed.set_values(0, std::move(values)); + + specmicp::Matrix numatrix(1, thedatabase->nb_component()); + numatrix << 0, 0, 2, 1; + thedatabase->sorbed.set_nu_matrix(numatrix); + + + Vector total_concentration = Vector::Zero(thedatabase->nb_component()); + total_concentration(id_h2o) = 55.5; + total_concentration(id_oh) = 2e-3; + total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.enable_conservation_water(); constraints.enable_surface_model(1.23456); specmicp::AdimensionalSystem system (thedatabase, constraints); - specmicp::Vector x(4); - x << 1.0, -3, -3, -0.2; + Vector x = Vector::Zero(system.total_dofs()); + x(system.ideq_w()) = 1.0; + x(system.ideq_paq(id_oh)) = std::log10(2e-3); + x(system.ideq_paq(id_ca)) = -3; + specmicp::Vector residual; system.get_residuals(x, residual); + std::cout << residual << std::endl; specmicp::Matrix analytical_jacobian; system.analytical_jacobian(x, analytical_jacobian); specmicp::Matrix fd_jacobian; system.finite_difference_jacobian(x, fd_jacobian); std::cout << analytical_jacobian - fd_jacobian << std::endl; REQUIRE((analytical_jacobian - fd_jacobian).isMuchSmallerThan(1.0, 1e-3)); } } diff --git a/tests/specmicp/adim/adimensional_system_carboalu.cpp b/tests/specmicp/adim/adimensional_system_carboalu.cpp index 59b4c2b..16b9f09 100644 --- a/tests/specmicp/adim/adimensional_system_carboalu.cpp +++ b/tests/specmicp/adim/adimensional_system_carboalu.cpp @@ -1,137 +1,144 @@ #include #include "catch.hpp" #include "utils/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 "utils/timer.hpp" + +#include "database/database.hpp" + TEST_CASE("Carboalu - using adimensional system ", "[Adimensional, Carbonation, carboalu]") { + std::cerr.flush(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Carboalu") { + std::cout << "---------------------\n Carboalu\n ---------------------" << std::endl; + + specmicp::Timer timer_carboalu; + timer_carboalu.start(); + specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"}, }); thedatabase.swap_components(swapping); + //thedatabase.remove_half_cell_reactions(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; specmicp::scalar_t mult = 5e3; specmicp::scalar_t m_c3s = mult*0.6; specmicp::scalar_t m_c2s = mult*0.2; specmicp::scalar_t m_c3a = mult*0.10; specmicp::scalar_t m_gypsum = mult*0.10; specmicp::scalar_t wc = 0.8; specmicp::scalar_t m_water = wc*1e-3*( m_c3s*(3*56.08+60.08) + m_c2s*(2*56.06+60.08) + m_c3a*(3*56.08+101.96) + m_gypsum*(56.08+80.06+2*18.02) ); formulation.mass_solution = m_water; formulation.amount_minerals = { {"C3S", m_c3s}, {"C2S", m_c2s}, {"C3A", m_c3a}, {"Gypsum", m_gypsum} }; formulation.extra_components_to_keep = {"HCO3[-]", }; formulation.minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Calcite", "Al(OH)3,am", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", "Straetlingite", "Gypsum", "Ettringite", "Thaumasite" }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); - specmicp::index_t id_h2o = 0; - specmicp::index_t id_ho = 1; - specmicp::index_t id_co2 = 3; + specmicp::index_t id_h2o = raw_data->get_id_component("H2O"); + specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); + specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; - options.solver_options.max_iter = 75; + options.solver_options.max_iter = 100; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-6; options.solver_options.steptol = 1e-10; options.system_options.non_ideality_tolerance = 1e-12; - specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral); + specmicp::Vector x; specmicp::scalar_t dh2co3 = mult*0.01; - for (specmicp::index_t i: raw_data->range_component()) - { - std::cout << raw_data->labels_basis[i] << std::endl; - } - specmicp::uindex_t tot_iter = 0; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.5, {}); - x.segment(1, raw_data->nb_component-1).setConstant(-6.0); - + solver.initialise_variables(x, 0.5, -3.0); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); - std::cout << perf.nb_iterations << std::endl; + std::cout << "First point number of iterations : " << perf.nb_iterations << " (Ref: 44)"<< std::endl; tot_iter += perf.nb_iterations; specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); for (int k=0; k<250; ++k) { constraints.total_concentrations(id_h2o) += dh2co3; constraints.total_concentrations(id_ho) -= dh2co3; constraints.total_concentrations(id_co2) += dh2co3; solver = specmicp::AdimensionalSystemSolver (raw_data, constraints, solution, options); //std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl; specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); //std::cout << x << std::endl; REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); //std::cout << perf.nb_iterations << std::endl; tot_iter += perf.nb_iterations; solution = solver.get_raw_solution(x); //specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); //std::cout << solution.main_variables(0) << std::endl; //std::cout << 14+solution.main_variables(1) << std::endl; } - std::cout << tot_iter << std::endl; + std::cout << "Number of iterations : " << tot_iter << " (Ref: 466)"<< std::endl; + timer_carboalu.stop(); + + std::cout << "Elapsed time : " << timer_carboalu.elapsed_time() << "s (Ref: 0.028s)" << std::endl; } } diff --git a/tests/specmicp/adim/adimensional_system_conditions.cpp b/tests/specmicp/adim/adimensional_system_conditions.cpp index ce11352..e2d9fc8 100644 --- a/tests/specmicp/adim/adimensional_system_conditions.cpp +++ b/tests/specmicp/adim/adimensional_system_conditions.cpp @@ -1,243 +1,248 @@ #include #include "catch.hpp" #include "utils/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/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" +#include "database/database.hpp" + TEST_CASE("Fixed activity", "[specmicp, MiCP, program, adimensional, solver, conditions]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; SECTION("Thermocarbo - saturated ") { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; specmicp::scalar_t mult = 1.0; formulation.mass_solution = mult*0.156; formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; constraints.water_equation = specmicp::WaterEquationType::SaturatedSystem; - specmicp::index_t id_h2o = 0; - specmicp::index_t id_ho = 1; - specmicp::index_t id_co2 = 2; + specmicp::index_t id_h2o = raw_data->water_index(); + specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); + specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; - options.solver_options.use_scaling = false; + options.solver_options.use_scaling = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 200; options.solver_options.condition_limit = -1; + options.solver_options.set_tolerance(1e-10,1e-12); + options.system_options.cutoff_total_concentration = 1e-16; - specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral+1); - x.setZero(); - x(0) = 0.8; - x.segment(1, raw_data->nb_component-1).setConstant(-2.0); + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); + + specmicp::Vector x; + solver.initialise_variables(x, 0.8, 2); specmicp::scalar_t dh2co3 = 0.1; for (int k=0; k<40; ++k) { specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); - constraints.total_concentrations(id_h2o) += mult*dh2co3; - constraints.total_concentrations(id_ho) -= mult*dh2co3; - constraints.total_concentrations(id_co2) += mult*dh2co3; specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set); //std::cout << solution.main_variables(0) << std::endl; - std::cout << 14+solution.main_variables(1) << std::endl; + std::cout << sol.pH() << std::endl; + constraints.total_concentrations(id_h2o) += mult*dh2co3; + constraints.total_concentrations(id_ho) -= mult*dh2co3; + constraints.total_concentrations(id_co2) += mult*dh2co3; } } SECTION("Carbonate speciation - fixed activity") { specmicp::database::Database dbhandler("../data/cemdata_specmicp.js"); std::vector to_keep = {"H2O", "H[+]", "HCO3[-]", "Na[+]"}; dbhandler.keep_only_components(to_keep); specmicp::RawDatabasePtr thedatabase = dbhandler.get_database(); specmicp::index_t id_h = dbhandler.component_label_to_id("H[+]"); specmicp::index_t id_na = dbhandler.component_label_to_id("Na[+]"); specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]"); specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2"); specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]"); specmicp::index_t id_naco3 = dbhandler.aqueous_label_to_id("NaCO3[-]"); specmicp::index_t id_nahco3 = dbhandler.aqueous_label_to_id("NaHCO3"); specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)"); - specmicp::Vector total_concentrations(thedatabase->nb_component); + specmicp::Vector total_concentrations(thedatabase->nb_component()); total_concentrations(0) = 55; total_concentrations(id_hco3) = 0.1; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = id_na; constraints.enable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.system_options.non_ideality = true; options.units_set.length = specmicp::units::LengthUnit::decimeter; options.solver_options.max_iter = 50; options.solver_options.maxstep = 10; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.non_monotone_linesearch = false; options.solver_options.penalization_factor = 0.8; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 10; options.solver_options.fvectol = 1e-6; specmicp::Vector x; constraints.add_fixed_activity_component(id_h, -4.0); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints, options); solver.initialise_variables(x, 0.8, -2); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << 4.0 << " \t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << " \t" << sol.molality_aqueous(id_naco3) << " \t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.gas_fugacity(id_co2g) << std::endl; for (double ph=4.5; ph<13; ph+=0.5) { constraints.fixed_activity_cs[0].log_value = -ph; solver = specmicp::AdimensionalSystemSolver( thedatabase, constraints, options); perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << ph << "\t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << "\t" << sol.molality_aqueous(id_naco3) << "\t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.gas_fugacity(id_co2g) << std::endl; } //std::cout << x << std::endl; } SECTION("Carbonate speciation - fixed fugacity") { specmicp::database::Database dbhandler("../data/cemdata_specmicp.js"); std::vector to_keep = {"H2O", "H[+]", "HCO3[-]"}; dbhandler.keep_only_components(to_keep); specmicp::RawDatabasePtr thedatabase = dbhandler.get_database(); - specmicp::index_t id_h = 1; - specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]"); + specmicp::index_t id_h = thedatabase->get_id_component("H[+]"); + specmicp::index_t id_hco3 = thedatabase->get_id_component("HCO3[-]"); - specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2"); - specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]");; - specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)"); + specmicp::index_t id_co2 = thedatabase->get_id_aqueous("CO2"); + specmicp::index_t id_co3 = thedatabase->get_id_aqueous("CO3[2-]");; + specmicp::index_t id_co2g = thedatabase->get_id_gas("CO2(g)"); - specmicp::Vector total_concentrations(thedatabase->nb_component); + specmicp::Vector total_concentrations(thedatabase->nb_component()); total_concentrations(0) = 55; total_concentrations(id_hco3) = 0.1; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = id_h; constraints.enable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.system_options.non_ideality = true; options.units_set.length = specmicp::units::LengthUnit::decimeter; options.solver_options.max_iter = 100; options.solver_options.maxstep = 10; options.solver_options.maxiter_maxstep = 10; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.penalization_factor = 0.8; options.solver_options.factor_descent_condition = -1; options.solver_options.fvectol = 1e-8; options.solver_options.condition_limit = -1; specmicp::Vector x; constraints.add_fixed_fugacity_gas(id_co2g, id_hco3, -5); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints, options); solver.initialise_variables(x, 0.8, -2.0); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << -5 << " \t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << " \t" << sol.gas_fugacity(id_co2g) << std::endl; for (double lfug=-4.5; lfug<=-0.5; lfug+=0.25) { constraints.fixed_fugacity_cs[0].log_value = lfug; solver = specmicp::AdimensionalSystemSolver( thedatabase, constraints, options); perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << lfug << "\t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << "\t" << sol.gas_fugacity(id_co2g) << std::endl; } //std::cout << x << std::endl; } } diff --git a/tests/specmicp/adim/adimensional_system_pcfm.cpp b/tests/specmicp/adim/adimensional_system_pcfm.cpp index 33d7e8f..be1d1cf 100644 --- a/tests/specmicp/adim/adimensional_system_pcfm.cpp +++ b/tests/specmicp/adim/adimensional_system_pcfm.cpp @@ -1,54 +1,61 @@ #include "catch.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp/adimensional/adimensional_system_pcfm.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" +#include "database/database.hpp" + specmicp::RawDatabasePtr get_pcfm_test_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_solid_phases(); return thedatabase.get_database(); } // Test the Positive continuous fraction method TEST_CASE("Positive continuous fraction method", "[specmicp, MiCP, program, adimensional, PCFM]") { specmicp::RawDatabasePtr the_database = get_pcfm_test_database(); + auto id_h2o = the_database->water_index(); + auto id_oh = the_database->get_id_component("HO[-]"); + auto id_ca = the_database->get_id_component("Ca[2+]"); SECTION("PCFM") { - specmicp::Vector total_concentration(the_database->nb_component); - total_concentration << 55.5, 2e-2, 1e-3; + specmicp::Vector total_concentration = specmicp::Vector::Zero(the_database->nb_component()); + total_concentration(id_h2o) = 55.5; + total_concentration(id_oh) = 2e-3; + total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.disable_conservation_water(); constraints.enable_surface_model(1.23456); - constraints.total_concentrations(2) = 0; + //constraints.total_concentrations(2) = 0; std::shared_ptr ptrsystem = std::make_shared(the_database, constraints); specmicp::AdimensionalSystemSolver solver(the_database, constraints); specmicp::AdimensionalSystemPCFM pcfm_solver(ptrsystem); specmicp::Vector x; solver.initialise_variables(x, 1.0, -3.0, 0.0); specmicp::PCFMReturnCode retcode = pcfm_solver.solve(x); REQUIRE(retcode == specmicp::PCFMReturnCode::Success); } } diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index 8c6732c..b2c246d 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,101 +1,104 @@ #include "catch.hpp" #include #include "utils/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/adimensional/adimensional_system_solution_extractor.hpp" +#include "database/database.hpp" 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("../data/cemdata_specmicp.js"); std::map 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 = 1; + 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 = 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+1); - x.setConstant(0.0); - x(0) = 0.8; - x.segment(1, raw_data->nb_component-1).setConstant(-2); + specmicp::Vector x; + solver.initialise_variables(x, 0.8, -2); - specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true); + specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); - REQUIRE(x(0) == Approx(0.802367).epsilon(1e-4)); - REQUIRE(x(4) == Approx(0.32947).epsilon(1e-2)); + 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("../data/cemdata_specmicp.js"); std::map 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 = {{"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 = 1; + 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 = false; 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::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); } } diff --git a/tests/specmicp/adim/adimensional_system_solver.cpp b/tests/specmicp/adim/adimensional_system_solver.cpp index 91b7da0..231ca2b 100644 --- a/tests/specmicp/adim/adimensional_system_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_solver.cpp @@ -1,148 +1,160 @@ #include "catch.hpp" #include #include "utils/log.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "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 "database/database.hpp" specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); - std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]"}; + std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); 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(); - SECTION("Solving simple case") { - - + 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"); - specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); - - specmicp::Vector total_concentrations(thedatabase->nb_component); - total_concentrations << 3.0e4, 2e4 , 1e4; + SECTION("Solving simple case") { - specmicp::AdimensionalSystemConstraints constraints(total_concentrations); + 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.use_crashing = false; - options.use_scaling = true; - options.penalization_factor = 0.8; - options.max_factorization_step = 4; - options.factor_descent_condition = 1e-10; + options.disable_crashing(); + options.disable_scaling(); + options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); - specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral); - x << 0.8, -2.0, -2.3, 0.0; + 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); - REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4)); - REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); - REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); - REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4)); + CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); + CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.48059).epsilon(1e-4)); + CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.8538).epsilon(1e-4)); + CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volum in cm^3") { - - specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); - - specmicp::Vector total_concentrations(thedatabase->nb_component); - total_concentrations << 0.03, 0.02, 0.01; - specmicp::AdimensionalSystemConstraints constraints(total_concentrations); + 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 = 100; + options.maxstep = 10; options.maxiter_maxstep = 100; - options.use_crashing = false; - options.use_scaling = false; - options.penalization_factor = 0.8; - options.max_factorization_step = 4; - options.factor_descent_condition = 1e-10; + options.disable_crashing(); + options.enable_scaling(); + options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); system->length_unit() = specmicp::units::LengthUnit::centimeter; specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); - specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral); - x << 0.8, -2.0, -2.3, 0.0; + 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); - REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4)); - REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); - REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); - REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4)); + CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); + CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.48059).epsilon(1e-4)); + CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.8538).epsilon(1e-4)); + CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Automatic solver") { - specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); - specmicp::Vector total_concentrations(thedatabase->nb_component); - total_concentrations << 0.03, 0.02, 0.01; - specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral+1); - x << 0.8, -2.0, -2.3, 0.0, 0.0; - specmicp::AdimensionalSystemConstraints constraints(total_concentrations); + 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.initialise_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 = 100.0; + 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 = false; - solver.get_options().solver_options.factor_descent_condition = 1e-10; + 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; solver.solve(x); + AdimensionalSystemSolution solution = solver.get_raw_solution(x); + AdimensionalSystemSolutionExtractor extr(solution, thedatabase, solver.get_options().units_set); - REQUIRE(x(0) == Approx(0.542049).epsilon(1e-4)); - REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); - REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); - REQUIRE(x(3) == -HUGE_VAL); - REQUIRE(x(4) == Approx(0.32966).epsilon(1e-4)); + CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4)); + CHECK(extr.log_molality_component(id_oh) == Approx(-1.48059).epsilon(1e-4)); + CHECK(extr.log_molality_component(id_ca) == Approx(-1.8538).epsilon(1e-4)); + //CHECK(extr.free_surface_concentration() == -HUGE_VAL); + CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4)); } } diff --git a/tests/specmicp/adim/adimensional_system_thermocarbo.cpp b/tests/specmicp/adim/adimensional_system_thermocarbo.cpp index b0120e2..03045ac 100644 --- a/tests/specmicp/adim/adimensional_system_thermocarbo.cpp +++ b/tests/specmicp/adim/adimensional_system_thermocarbo.cpp @@ -1,72 +1,144 @@ #include "catch.hpp" #include #include "utils/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 "utils/timer.hpp" + +#include "database/database.hpp" + TEST_CASE("thermocarbo - using adimensional system ", "[Adimensional, Thermocarbo]") { + std::cerr.flush(); + std::cout.flush(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Thermocarbo") { + specmicp::Timer tot_timer; + tot_timer.start(); + + std::cout << "\n-------------------\n Thermocarbo\n -------------------" << std::endl; + + specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); + + specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); + std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); - specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; specmicp::scalar_t mult = 1.0; formulation.mass_solution = mult*0.156; formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); - constraints.charge_keeper = 1; + constraints.charge_keeper = raw_data->get_id_component("HO[-]"); - specmicp::index_t id_h2o = 0; - specmicp::index_t id_ho = 1; - specmicp::index_t id_co2 = 2; + specmicp::index_t id_h2o = raw_data->water_index(); + specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); + specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 10.0; options.solver_options.maxiter_maxstep = 100; - options.solver_options.use_crashing = false; - options.solver_options.use_scaling = false; - options.solver_options.factor_descent_condition = -1; - options.solver_options.factor_gradient_search_direction = 200; - - specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral+1); - x.setZero(); - x(0) = 0.8; - x.segment(1, raw_data->nb_component-1).setConstant(-2.0); + options.solver_options.disable_descent_direction(); + options.solver_options.set_tolerance(1e-8, 1e-12); + options.solver_options.disable_non_monotone_linesearch(); + options.solver_options.disable_condition_check(); + options.solver_options.disable_crashing(); + options.solver_options.disable_scaling(); + + + specmicp::Vector x; + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); + solver.initialise_variables(x, 0.8, -4.0); specmicp::scalar_t dh2co3 = 0.1; - for (int k=0; k<40; ++k) + + specmicp::index_t nb_iter {0}; + + constexpr specmicp::index_t nb_step = 40; + specmicp::Vector solution(nb_step); + + for (int k=0; k= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); + REQUIRE(static_cast(perf.return_code) >= + static_cast(specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)); + nb_iter += perf.nb_iterations; constraints.total_concentrations(id_h2o) += mult*dh2co3; constraints.total_concentrations(id_ho) -= mult*dh2co3; constraints.total_concentrations(id_co2) += mult*dh2co3; - specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); - //std::cout << solution.main_variables(0) << std::endl; - std::cout << 14+solution.main_variables(1) << std::endl; + specmicp::AdimensionalSystemSolution adim_solution = solver.get_raw_solution(x); + solution(k) = 14+adim_solution.main_variables(id_ho); + } + + tot_timer.stop(); + + std::cout << "Total time : " << tot_timer.elapsed_time() << "s (Ref 0.040s)" << std::endl; + std::cout << "Nb iter : " << nb_iter << " (Ref 1846)" << std::endl; + specmicp::Vector reference_solution(nb_step); + reference_solution << 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.5194, + 12.1374, + 12.1374, + 12.1374, + 12.1374, + 12.1374, + 12.1374, + 12.1374, + 12.1374, + 9.8575, + 9.8575, + 9.8575, + 9.8575, + 9.8575, + 9.8575, + 9.8575, + 9.8575, + 8.97554, + 5.33701, + 5.15035, + 5.0435, + 4.96886, + 4.91174, + 4.86562, + 4.82705, + 4.794, + 4.76515, + 4.73958, + 4.71668, + 4.69597; + for (int k=0;k