diff --git a/src/specmicp/problem_solver/reactant_box.cpp b/src/specmicp/problem_solver/reactant_box.cpp index 1f0c3ba..155f22d 100644 --- a/src/specmicp/problem_solver/reactant_box.cpp +++ b/src/specmicp/problem_solver/reactant_box.cpp @@ -1,563 +1,572 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "reactant_box.hpp" #include "specmicp/adimensional/adimensional_system_structs.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp_database/database.hpp" #include "specmicp_database/unknown_class.hpp" #include "specmicp_common/physics/laws.hpp" #include #include #include #include namespace specmicp { struct SPECMICP_DLL_LOCAL fug_container { std::string gas; std::string component; scalar_t fugacity; fug_container( const std::string& agas, const std::string& acomponent, scalar_t t_fugacity ): gas(agas), component(acomponent), fugacity(t_fugacity) {} }; struct ReactantBox::ReactantBoxImpl { scalar_t mass_solution {-1}; std::vector keep_components {}; std::unordered_map aqueous_species {}; std::unordered_map solid_phase {}; // constraints std::string charge_keeper {""}; scalar_t inert_vol_frac {0.0}; WaterEquationType water_eq_type {WaterEquationType::MassConservation}; water_partial_pressure_f h2o_partial_pressure_model; std::vector fix_fugacity; std::vector> fix_activity; std::vector> fix_molality; std::vector get_components_to_keep(std::shared_ptr) const; void dissolve_aqueous_species( std::shared_ptr raw_data, Vector& total_concentration ) const; void dissolve_solid_phases( std::shared_ptr raw_data, Vector& total_concentration ) const; void raise_unknown_unit(const std::string& error_msg) { throw std::invalid_argument(error_msg); } void assert_mass_solution() { if (mass_solution == -1) { throw std::logic_error("Mass of solution is not set !"); } } scalar_t factor_total_concentration(const units::UnitsSet& units_set) const; }; ReactantBox::ReactantBox( std::shared_ptr raw_data, const units::UnitsSet& units_set ): units::UnitBaseClass(units_set), database::DatabaseHolder(raw_data), m_impl(utils::make_pimpl()) { } ReactantBox::ReactantBox(ReactantBox&& other): units::UnitBaseClass(other.get_units()), database::DatabaseHolder(other.m_data), m_impl(std::move(other.m_impl)) { } ReactantBox& ReactantBox::operator=(ReactantBox&& other) { m_impl = std::move(other.m_impl); return *this; } ReactantBox::~ReactantBox() = default; void ReactantBox::set_solution(scalar_t value, std::string str_unit) { set_solution(value, units::parse_amount_unit(str_unit)); } void ReactantBox::set_solution( scalar_t value, const units::AmountUnit& unit ) { scalar_t factor; if (unit.type == units::AmountUnitType::Mass) { factor = unit.factor_si; } else if (unit.type == units::AmountUnitType::Volume) { factor = unit.factor_si*laws::density_water(units::SI_units); + } else if (unit.type == units::AmountUnitType::VolumeFraction) { + factor = unit.factor_si*laws::density_water(units::SI_units); } else if (unit.type == units::AmountUnitType::MassConcentration) { factor = unit.factor_si; } else { m_impl->raise_unknown_unit( - "Only mass. volume, or mass concentration units " + "Only mass. volume, volume fraction," + " or mass concentration units " "are accepted for the amount of solution" ); } m_impl->mass_solution = factor*value; } void ReactantBox::add_aqueous_species( std::string name, scalar_t value, std::string str_unit ) { add_aqueous_species(name, value, units::parse_amount_unit(str_unit)); } void ReactantBox::add_aqueous_species( std::string name, scalar_t value, const units::AmountUnit& unit ) { scalar_t factor = 1.0; if (unit.type == units::AmountUnitType::Molality) { factor = unit.factor_si; } else if (unit.type == units::AmountUnitType::MoleConcentration) { factor = unit.factor_si/laws::density_water(units::SI_units); } else if (unit.type == units::AmountUnitType::MassConcentration) { const scalar_t molar_mass = database::molar_mass_aqueous(get_database(), name); factor = unit.factor_si / molar_mass / laws::density_water(units::SI_units); } else if (unit.type == units::AmountUnitType::Mass) { m_impl->assert_mass_solution(); const scalar_t molar_mass = database::molar_mass_aqueous(get_database(), name); factor = unit.factor_si / molar_mass / m_impl->mass_solution; } else if (unit.type == units::AmountUnitType::NumberOfMoles) { m_impl->assert_mass_solution(); factor = unit.factor_si / m_impl->mass_solution; } else { m_impl->raise_unknown_unit("Unit for aqueous species " + name + "is invalid." " Only accepted units for aqueous species are" " molality, mole and mass concettration," " mass and number of moles."); } const scalar_t molality = factor*value; m_impl->aqueous_species[name] = molality; } void ReactantBox::add_solid_phase( std::string name, scalar_t value, std::string str_unit ) { return add_solid_phase(name, value, units::parse_amount_unit(str_unit)); } void ReactantBox::add_solid_phase( std::string name, scalar_t value, const units::AmountUnit& unit ) { scalar_t factor = 1.0; // convert to mole concentration is SI if (unit.type == units::AmountUnitType::MoleConcentration or unit.type == units::AmountUnitType::NumberOfMoles) { factor = unit.factor_si; } else if (unit.type == units::AmountUnitType::MassConcentration or unit.type == units::AmountUnitType::Mass ) { const scalar_t molar_mass = database::molar_mass_solid_phase(get_database(), name); factor = unit.factor_si / molar_mass; } else if (unit.type == units::AmountUnitType::Volume) { const scalar_t molar_volume = database::molar_volume_solid_phase(get_database(), name); factor = unit.factor_si / molar_volume; } + else if (unit.type == units::AmountUnitType::VolumeFraction) + { + const scalar_t molar_volume = database::molar_volume_solid_phase(get_database(), name); + factor = unit.factor_si / molar_volume; + + } else { m_impl->raise_unknown_unit("Unit for solid phase " + name + " is invalid."); } m_impl->solid_phase[name] = factor*value; } void ReactantBox::add_component(const std::string& name) { m_impl->keep_components.push_back(name); } Vector ReactantBox::get_total_concentration(bool modify_db) const { std::shared_ptr raw_data = get_database(); if (modify_db) { // remove extra components const auto id_comp_to_keep = m_impl->get_components_to_keep(raw_data); database::Database db_manager(raw_data); db_manager.keep_only_components(id_comp_to_keep); } Vector total_concentration = Vector::Zero(raw_data->nb_component()); total_concentration[database::water_id] = m_impl->mass_solution/raw_data->molar_mass_basis(database::water_id, units::SI_units); m_impl->dissolve_aqueous_species(raw_data, total_concentration); m_impl->dissolve_solid_phases(raw_data, total_concentration); total_concentration *= m_impl->factor_total_concentration(get_units()); return total_concentration; } void ReactantBox::ReactantBoxImpl::dissolve_aqueous_species( std::shared_ptr raw_data, Vector& total_concentration ) const { for (auto& it: aqueous_species) { database::GenericAqueousSpecies aq = database::find_aqueous_species(raw_data, it.first); const scalar_t conc = mass_solution*it.second; if (aq.type == database::AqueousSpeciesClass::Component) { total_concentration[aq.id] += conc; } else if (aq.type == database::AqueousSpeciesClass::Aqueous) { for (auto idc: raw_data->range_component()) { const auto nu = raw_data->nu_aqueous(aq.id, idc); if (nu != 0) { total_concentration[idc] += nu*conc; } } } else if (aq.type == database::AqueousSpeciesClass::Compound) { for (auto idc: raw_data->range_component()) { const auto nu = raw_data->nu_compound(aq.id, idc); if (nu != 0) { total_concentration[idc] += nu*conc; } } } else { throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species."); } } } void ReactantBox::ReactantBoxImpl::dissolve_solid_phases( std::shared_ptr raw_data, Vector& total_concentration ) const { for (auto& it: solid_phase) { database::GenericSolidPhase sp = database::find_solid_phase(raw_data, it.first); if (sp.type == database::SolidPhaseClass::EquilibriumMineral) { for (auto idc: raw_data->range_component()) { const auto nu = raw_data->nu_mineral(sp.id, idc); if (nu != 0) { total_concentration[idc] += nu*it.second; } } } else if (sp.type == database::SolidPhaseClass::MineralKinetics) { for (auto idc: raw_data->range_component()) { const auto nu = raw_data->nu_mineral_kinetic(sp.id, idc); if (nu != 0) { total_concentration[idc] += nu*it.second; } } } else { throw std::invalid_argument("'"+it.first+"' is not a valid solid phase."); } } } std::vector ReactantBox::ReactantBoxImpl::get_components_to_keep( std::shared_ptr raw_data ) const { std::vector to_keep(raw_data->nb_component(), false); to_keep[database::water_id] = true; to_keep[database::electron_id] = true; // aqueous species // --------------- for (auto& it: aqueous_species) { database::GenericAqueousSpecies aq = database::find_aqueous_species(raw_data, it.first); if (aq.type == database::AqueousSpeciesClass::Component) { to_keep[aq.id] = true; } else if (aq.type == database::AqueousSpeciesClass::Aqueous) { for (auto idc: raw_data->range_aqueous_component()) { if (raw_data->nu_aqueous(aq.id, idc) != 0) { to_keep[idc] = true; } } } else if (aq.type == database::AqueousSpeciesClass::Compound) { for (auto idc: raw_data->range_aqueous_component()) { if (raw_data->nu_compound(aq.id, idc) != 0) { to_keep[idc] = true; } } } else { throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species."); } } // solid phases // ------------ for (auto& it: solid_phase) { database::GenericSolidPhase sp = database::find_solid_phase(raw_data, it.first); if (sp.type == database::SolidPhaseClass::EquilibriumMineral) { for (auto idc: raw_data->range_aqueous_component()) { if (raw_data->nu_mineral(sp.id, idc) != 0) { to_keep[idc] = true; } } } else if (sp.type == database::SolidPhaseClass::MineralKinetics) { for (auto idc: raw_data->range_aqueous_component()) { if (raw_data->nu_mineral_kinetic(sp.id, idc) != 0) { to_keep[idc] = true; } } } else { throw std::invalid_argument("'"+it.first+"' is not a valid solid phase."); } } // extra components // ----------------- for (auto& comp_label: keep_components) { const index_t idc = raw_data->get_id_component(comp_label); to_keep[idc] = true; } // get ids of components to keep std::vector id_comp_to_keep; id_comp_to_keep.reserve(raw_data->nb_component()); for (auto id: raw_data->range_component()) { if (to_keep[id]) { id_comp_to_keep.push_back(id); } } return id_comp_to_keep; } scalar_t ReactantBox::ReactantBoxImpl::factor_total_concentration( const units::UnitsSet& units_set ) const { scalar_t mol_fac = units::scaling_factor(units_set.quantity); scalar_t length_fac = units::scaling_factor(units_set.length); return std::pow(length_fac, 3)/mol_fac; } void ReactantBox::set_charge_keeper(const std::string& charge_keeper) { m_impl->charge_keeper = charge_keeper; } void ReactantBox::set_inert_volume_fraction(scalar_t inert_volume_fraction) { m_impl->inert_vol_frac = inert_volume_fraction; } void ReactantBox::set_saturated_system() { m_impl->water_eq_type = WaterEquationType::SaturatedSystem; } void ReactantBox::disable_conservation_water() { m_impl->water_eq_type = WaterEquationType::NoEquation; } void ReactantBox::set_water_partial_pressure_model( water_partial_pressure_f h2o_pressure_model ) { m_impl->h2o_partial_pressure_model = h2o_pressure_model; } void ReactantBox::add_fixed_fugacity_gas( std::string gas, std::string component, scalar_t fugacity ) { m_impl->fix_fugacity.emplace_back(gas, component, fugacity); add_component(component); } void ReactantBox::add_fixed_activity_component( std::string component, scalar_t activity ) { m_impl->fix_activity.emplace_back(component, activity); add_component(component); } void ReactantBox::add_fixed_molality_component( std::string component, scalar_t molality ) { m_impl->fix_molality.emplace_back(component, molality); add_component(component); } AdimensionalSystemConstraints ReactantBox::get_constraints(bool modify_db) const { database::DataContainer* raw_data = get_database().get(); AdimensionalSystemConstraints constraints; constraints.total_concentrations = get_total_concentration(modify_db); if (not m_impl->charge_keeper.empty()) { constraints.set_charge_keeper( raw_data->get_id_component(m_impl->charge_keeper) ); } constraints.set_inert_volume_fraction( m_impl->inert_vol_frac ); constraints.water_equation = m_impl->water_eq_type; if (m_impl->h2o_partial_pressure_model != nullptr) { constraints.set_water_partial_pressure_model(m_impl->h2o_partial_pressure_model); } for (auto& fug: m_impl->fix_fugacity) { constraints.add_fixed_fugacity_gas( raw_data->get_id_gas(fug.gas), raw_data->get_id_component(fug.component), std::log10(fug.fugacity) ); } for (auto& act: m_impl->fix_activity) { constraints.add_fixed_activity_component( raw_data->get_id_component(act.first), std::log10(act.second) ); } for (auto& mol: m_impl->fix_molality) { constraints.add_fixed_molality_component( raw_data->get_id_component(mol.first), std::log10(mol.second) ); } return constraints; } } // end namespace specmicp diff --git a/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp index ea9daba..0e9ac23 100644 --- a/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp +++ b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp @@ -1,152 +1,155 @@ #include "catch.hpp" #include "reactmicp/io/hdf5_unsaturated.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp_database/database.hpp" #include "specmicp/problem_solver/reactant_box.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "reactmicp/systems/unsaturated/variables_interface.hpp" #include "dfpm/mesh.hpp" #include #define PATH_TEST_FILES "test_hdf5_unsaturated" using namespace specmicp; using namespace specmicp::reactmicp::systems::unsaturated; static specmicp::database::RawDatabasePtr get_database() { static database::RawDatabasePtr raw_data {nullptr}; if (raw_data == nullptr) { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); thedatabase.keep_only_components( {"H2O", "H[+]", "Ca[2+]", "Si(OH)4", "HCO3[-]"} ); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } static AdimensionalSystemSolution get_init_solution( database::RawDatabasePtr the_database, units::UnitsSet the_units ) { specmicp::ReactantBox reactant_box(the_database, the_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_solution(0.5, "L/L"); + reactant_box.add_solid_phase("C3S", 0.350, "L/L"); + reactant_box.add_solid_phase("C2S", 0.50, "L/L"); reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); reactant_box.set_charge_keeper("HO[-]"); auto constraints = reactant_box.get_constraints(true); AdimensionalSystemSolver adim_solver(the_database, constraints); adim_solver.get_options().solver_options.set_tolerance(1e-14); adim_solver.get_options().units_set = the_units; Vector x; adim_solver.initialise_variables(x, 0.8, -3); micpsolver::MiCPPerformance perf = adim_solver.solve(x); - specmicp_assert(perf.return_code > micpsolver::MiCPSolverReturnCode::Success); + if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) + { + throw std::runtime_error("Error solving the initial problem"); + } return adim_solver.get_raw_solution(x); } TEST_CASE("hdf5 output", "[hdf5]") { std::string filename = "test_hdf5_unsaturated.hdf5"; database::RawDatabasePtr raw_data = get_database(); mesh::Uniform1DMeshGeometry mesh_geom; mesh_geom.nb_nodes = 10; mesh_geom.dx = 1; mesh_geom.section = 1; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(mesh_geom); units::UnitsSet the_units; the_units.length = units::LengthUnit::centimeter; specmicp::AdimensionalSystemSolution init_sol = get_init_solution(raw_data, the_units); specmicp::AdimensionalSystemSolutionExtractor extr(init_sol, raw_data, the_units); SECTION("Writer") { VariablesInterface vars_init(the_mesh, raw_data, {raw_data->get_id_component_from_element("C")}, the_units); for (auto node: the_mesh->range_nodes()) { vars_init.initialize_variables(node, extr); } auto vars = vars_init.get_variables(); io::UnsaturatedHDF5Saver saver(filename, vars, the_units); saver.save_timestep(10.0); } SECTION("Reader") { specmicp::io::UnsaturatedHDF5Reader reader(filename); // units units::UnitsSet read_unit = reader.get_units(); CHECK(read_unit.length == the_units.length); CHECK(read_unit.mass == the_units.mass); CHECK(read_unit.quantity == the_units.quantity); // adimensional solution specmicp::AdimensionalSystemSolution sol = reader.get_adim_solution(10.0, 1); REQUIRE(sol.is_valid); for (auto r: range(sol.main_variables.rows())) { if (std::isfinite(sol.main_variables(r))) { CHECK(sol.main_variables(r) == Approx(init_sol.main_variables(r)).epsilon(1e-10)); } } CHECK(sol.log_gamma.norm() == Approx(init_sol.log_gamma.norm()).epsilon(1e-10)); CHECK(sol.gas_fugacities.norm() == Approx(init_sol.gas_fugacities.norm()).epsilon(1e-10)); CHECK(sol.secondary_molalities.norm() == Approx(init_sol.secondary_molalities.norm()).epsilon(1e-10)); // main variables Vector liq_sat = reader.main_variable_vs_x(10.0, "H2O", "liquid_saturation"); CHECK(liq_sat.rows() == the_mesh->nb_nodes()); CHECK(liq_sat(1) == Approx(extr.saturation_water())); Vector ca_s = reader.main_variable_vs_x(10.0, "Ca[2+]", "solid_concentration"); CHECK(ca_s.rows() == the_mesh->nb_nodes()); CHECK(ca_s(1) == Approx(extr.total_solid_concentration(raw_data->get_id_component("Ca[2+]")))); } } diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index 45c8a8e..3839d0c 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,400 +1,467 @@ #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 +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 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 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 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 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 swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 15, "mol/kg"); reactant_box.set_charge_keeper("Cl[-]"); reactant_box.add_fixed_activity_component("H[+]", 1e-8); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.8, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.pH() == Approx(8)); } SECTION("Reactant box fixed molality") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 15, "mol/kg"); reactant_box.set_charge_keeper("Cl[-]"); reactant_box.add_fixed_molality_component("H[+]", 1e-8); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; solver.initialize_variables(x, 0.8, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.molality_component(raw_data->get_id_component("H[+]")) == Approx(1e-8)); } SECTION("Reactant box fixed fugacity") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map 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, -4); solver.get_options().solver_options.set_tolerance(1e-6, 1e-8); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.fugacity_gas(raw_data->get_id_gas("CO2(g)")) == Approx(1e-3)); } }