diff --git a/src/bin/specmicp/specmicp.cpp b/src/bin/specmicp/specmicp.cpp index 9450285..9629b5c 100644 --- a/src/bin/specmicp/specmicp.cpp +++ b/src/bin/specmicp/specmicp.cpp @@ -1,257 +1,277 @@ /* ============================================================================= 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 "specmicp_common/io/config_yaml_sections.h" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_database/database.hpp" #include "specmicp_database/io/configuration.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp/problem_solver/reactant_box.hpp" +#include "specmicp/problem_solver/smart_solver.hpp" + #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_saver.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/adimensional/adimensional_system_structs.hpp" #include "specmicp/io/configuration.hpp" #include "specmicp/io/print.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/io/configuration.hpp" #include #define SECTION_MAIN "__main__" #define SECTION_SPECIATION "speciation" #define SECTION_FORMULATION "formulation" #define SECTION_CONSTRAINTS "constraints" #define ATTRIBUTE_NAME "name" #define ATTRIBUTE_OUTPUT "output" #define ATTRIBUTE_OUTPUT_DB "output_db" namespace specmicp { void check_db(const std::string& filepath); void run_specmicp(io::YAMLConfigHandle& configuration); void run_formulation( io::YAMLConfigHandle&& conf_formulation, const RawDatabasePtr& the_database, const AdimensionalSystemSolverOptions& the_options ); } const static char* header = "SpecMiCP: Speciation solver based on the complementarity condition\n" "(c) Copyright 2014-2016 fabien Georget \n"; const static char* description = R"plop( Usage: specmicp specmicp -C/--check-db Options : the configuration file -c , --check-db check that the database can be parsed without error -h, --help produce this help )plop"; const static char* error_print_help = "Use -h to print help"; int main( int argc, const char* argv[]) { specmicp::init_logger(&std::cerr, specmicp::logger::Warning); if (argc < 2) { std::cerr << "Error: expecting one argument. "; std::cerr << error_print_help << std::endl; return EXIT_FAILURE; } std::cout << header << std::endl; const std::string first_arg = argv[1]; if (first_arg[0] == '-') { if (first_arg == "-h" or first_arg == "--help") { std::cout << description << std::endl; return EXIT_SUCCESS; } else if (first_arg == "-C" or first_arg == "--check-db") { if (argc < 3) { std::cerr << "Error : " << first_arg << " option requires an extra argument" << std::endl; return EXIT_FAILURE; } else { specmicp::check_db(argv[2]); return EXIT_SUCCESS; } } } const std::string conf_file = first_arg; auto conf = specmicp::io::YAMLConfigFile::load(conf_file); std::unique_ptr out_log; std::unique_ptr conf_log; if (conf.has_section(SPC_CF_S_LOGS)) { out_log = specmicp::io::configure_log(conf.get_section(SPC_CF_S_LOGS)); } if (conf.has_section(SPC_CF_S_CONF_LOGS)) { conf_log = specmicp::io::configure_conf_log(conf.get_section(SPC_CF_S_CONF_LOGS)); } specmicp::run_specmicp(conf); return EXIT_SUCCESS; } namespace specmicp { RawDatabasePtr get_db(io::YAMLConfigHandle& configuration) { return io::configure_database(configuration.get_section(SPC_CF_S_DATABASE)); } units::UnitsSet get_units(io::YAMLConfigHandle& configuration) { return io::configure_units(configuration.get_section(SPC_CF_S_UNITS)); } void run_specmicp(io::YAMLConfigHandle& configuration) { RawDatabasePtr the_database = get_db(configuration); units::UnitsSet the_units = get_units(configuration); AdimensionalSystemSolverOptions the_options; io::configure_specmicp_options(the_options, the_units, configuration.get_section(SPC_CF_S_SPECMICP)); if (not configuration.has_section(SECTION_SPECIATION)) { configuration.report_error(io::YAMLConfigError::MissingRequiredSection, "No system to solve !"); } auto conf_formulation = configuration.get_section(SECTION_SPECIATION); if (conf_formulation.is_sequence()) { uindex_t size = conf_formulation.size(); for (uindex_t ind=0; ind(ATTRIBUTE_NAME); std::string save_db = ""; if (conf.has_attribute(ATTRIBUTE_OUTPUT_DB)) { save_db = conf.get_attribute(ATTRIBUTE_OUTPUT_DB); } std::string save_solution = ""; if (conf.has_attribute(ATTRIBUTE_OUTPUT)) { save_solution = conf.get_attribute(ATTRIBUTE_OUTPUT); } else { WARNING << "No output file set, the solution will be lost in the void."; } auto rbox = io::configure_specmicp_reactant_box( the_database, the_options.units_set, std::move(conf) ); auto constraints = rbox.get_constraints(true); if (not save_db.empty()) { database::Database(the_database).save(save_db); } - AdimensionalSystemSolution solution = solve_equilibrium(the_database, constraints, the_options); + + SmartAdimSolver solver(the_database, constraints, the_options); + + // initialization + if (conf.has_section(SPC_CF_S_INITIALIZATION)) + { + io::configure_smart_solver_initialization( + solver, + conf.get_section(SPC_CF_S_INITIALIZATION) + ); + } + // solve the problen + + + auto is_solved = solver.solve(); + if (not is_solved) { + CRITICAL << "Failed to solve formulation : " << name << "."; + return; + } if (not save_solution.empty()) { AdimensionalSystemSolutionSaver(the_database).save_solution( - save_solution, solution); + save_solution, solver.get_solution()); } } void check_db(const std::string& filepath) { std::cout << "Checking database : " << filepath << "\n ------ \n"; database::Database the_database(filepath); if (not the_database.is_valid()) { throw std::invalid_argument("The database is not valid !"); } std::cout << "Database is valid !" << std::endl; } } //end namespace specmicp diff --git a/src/specmicp/io/configuration.cpp b/src/specmicp/io/configuration.cpp index ea3c717..3b59a74 100644 --- a/src/specmicp/io/configuration.cpp +++ b/src/specmicp/io/configuration.cpp @@ -1,742 +1,804 @@ /* ============================================================================= 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 "configuration.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/adimensional/config_default_options_solver.h" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp_common/physics/laws.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/reactant_box.hpp" +#include "specmicp/problem_solver/smart_solver.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/io/config_yaml_sections.h" +#include + namespace specmicp { namespace io { // Additional declarations // ======================= void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_activity, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_fugacity( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_fugacity, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_molality( AdimensionalSystemConstraints& constraints, const YAML::Node& conf_constraints_fixed_molality, const RawDatabasePtr& raw_db ); void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_activity ); void configure_specmicp_constraints_fixed_fugacity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_fugacity ); void configure_specmicp_constraints_fixed_molality( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_molality ); // Implementation // ============== void configure_specmicp_options( AdimensionalSystemSolverOptions& options, const units::UnitsSet& the_units, YAMLConfigHandle&& conf ) { micpsolver::MiCPSolverOptions& solver = options.solver_options; conf.set_if_attribute_exists( solver.fvectol, SPC_CF_S_SPECMICP_A_RES_TOL, 0.0); conf.set_if_attribute_exists( solver.steptol, SPC_CF_S_SPECMICP_A_STEP_TOL, 0.0); conf.set_if_attribute_exists( solver.max_iter, SPC_CF_S_SPECMICP_A_MAX_ITER, 0); conf.set_if_attribute_exists( solver.maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH, 0.0); conf.set_if_attribute_exists( solver.maxiter_maxstep, SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER, 0); conf.set_if_attribute_exists( solver.use_scaling, SPC_CF_S_SPECMICP_A_SCALING); conf.set_if_attribute_exists( solver.non_monotone_linesearch, SPC_CF_S_SPECMICP_A_NONMONOTONE); conf.set_if_attribute_exists( solver.factor_descent_condition, SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION); conf.set_if_attribute_exists( solver.condition_limit, SPC_CF_S_SPECMICP_A_COND_CHECK); conf.set_if_attribute_exists( solver.threshold_cycling_linesearch, SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH); AdimensionalSystemOptions& system = options.system_options; conf.set_if_attribute_exists( system.non_ideality_tolerance, SPC_CF_S_SPECMICP_A_NONIDEAL_TOL, 0.0); conf.set_if_attribute_exists( system.non_ideality_max_iter, SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER, 0); conf.set_if_attribute_exists( system.cutoff_total_concentration, SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC, 0.0); conf.set_if_attribute_exists( system.restart_concentration, SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION); conf.set_if_attribute_exists( system.restart_water_volume_fraction, SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC, 0.0); conf.set_if_attribute_exists( system.under_relaxation_factor, SPC_CF_S_SPECMICP_A_UNDER_RELAXATION, 0.0, 1.0); options.units_set = the_units; } // ReactantBox // ============ namespace internal { // helper functions static void configure_formulation( ReactantBox& reactant_box, YAMLConfigHandle&& configuration ); static void configure_formulation_aqueous( ReactantBox& reactant_box, YAMLConfigHandle&& cf_aqueous ); static void configure_formulation_solid( ReactantBox& reactant_box, YAMLConfigHandle&& cf_aqueous ); static void configure_constraints( ReactantBox& reactant_box, YAMLConfigHandle&& configuration ); static void configure_constraints_fixed_activity( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_act ); static void configure_constraints_fixed_fugacity( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_act ); static void configure_constraints_fixed_molality( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_mol ); } // end namespace internal ReactantBox SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box( RawDatabasePtr raw_db, const units::UnitsSet& the_units, YAMLConfigHandle&& configuration ) { ReactantBox reactant_box(raw_db, the_units); internal::configure_formulation( reactant_box, configuration.get_section(SPC_CF_S_FORMULATION) ); internal::configure_constraints( reactant_box, configuration.get_section(SPC_CF_S_CONSTRAINTS) ); return reactant_box; } namespace internal { // helper functions static void configure_formulation( ReactantBox& reactant_box, YAMLConfigHandle&& cf_formulation) { // solution auto cf_solution = cf_formulation.get_section(SPC_CF_S_FORMULATION_A_SOLUTION); const auto amount = cf_solution.get_required_attribute(SPC_CF_S_FORMULATION_A_AMOUNT); const auto unit = cf_solution.get_required_attribute(SPC_CF_S_FORMULATION_A_UNIT); reactant_box.set_solution(amount, unit); if (cf_formulation.has_section(SPC_CF_S_FORMULATION_A_AQUEOUS)) { configure_formulation_aqueous( reactant_box, cf_formulation.get_section(SPC_CF_S_FORMULATION_A_AQUEOUS) ); } if (cf_formulation.has_section(SPC_CF_S_FORMULATION_A_MINERALS)) { configure_formulation_solid( reactant_box, cf_formulation.get_section(SPC_CF_S_FORMULATION_A_MINERALS) ); } } static void configure_formulation_aqueous( ReactantBox& reactant_box, YAMLConfigHandle&& cf_aqueous ) { if (not cf_aqueous.is_sequence()) { cf_aqueous.report_error( YAMLConfigError::ListExpected, "Aqueous species must be provide as a list of triplet" " {label amount, unit}" ); } for (auto ind: RangeIterator(cf_aqueous.size())) { auto cf = cf_aqueous.get_section(ind); const auto label = cf.get_required_attribute(SPC_CF_S_FORMULATION_A_LABEL); const auto value = cf.get_required_attribute(SPC_CF_S_FORMULATION_A_AMOUNT); const auto unit = cf.get_required_attribute(SPC_CF_S_FORMULATION_A_UNIT); reactant_box.add_aqueous_species(label, value, unit); } } static void configure_formulation_solid( ReactantBox& reactant_box, YAMLConfigHandle&& cf_solid ) { if (not cf_solid.is_sequence()) { cf_solid.report_error( YAMLConfigError::ListExpected, "Solid species must be provide as a list of triplet" " {label amount, unit}" ); } for (auto ind: RangeIterator(cf_solid.size())) { auto cf = cf_solid.get_section(ind); const auto label = cf.get_required_attribute< std::string>(SPC_CF_S_FORMULATION_A_LABEL); const auto value = cf.get_required_attribute< double>(SPC_CF_S_FORMULATION_A_AMOUNT); const auto unit = cf.get_required_attribute< std::string>(SPC_CF_S_FORMULATION_A_UNIT); reactant_box.add_solid_phase(label, value, unit); } } static void configure_constraints( ReactantBox& reactant_box, YAMLConfigHandle&& cf) { // saturated system // ---------------- if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) { auto val = cf.get_attribute< bool>(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM); // only do something if the option is true if (val) { // check that conservation of water is not set auto val_conv = cf.get_optional_attribute( SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER, false); if (val_conv) { cf.report_error( YAMLConfigError::InvalidArgument, "saturated_system and Conservation_water are both" " set to true. Choose only one of them." ); } // everything is fine, we can set the option reactant_box.set_saturated_system(); } } // water conservation // ------------------ if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER)) { auto val = cf.get_attribute< bool>(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER); // we already checked that saturated_system and conservation water are // not both set if (not val) { // conservation of water is the default reactant_box.disable_conservation_water(); } } // Charge keeper // ------------- if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER)) { auto label = cf.get_attribute< std::string>(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER); reactant_box.set_charge_keeper(label); } // Fixed activity // --------------- if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY)) { configure_constraints_fixed_activity( reactant_box, cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY) ); } // Fixed fugacity // -------------- if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY)) { configure_constraints_fixed_fugacity( reactant_box, cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY) ); } // Fixed molality // --------------- if (cf.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY)) { configure_constraints_fixed_molality( reactant_box, cf.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY) ); } if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION)) { reactant_box.set_inert_volume_fraction( cf.get_attribute(SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION) ); } } static void configure_constraints_fixed_activity( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_act ) { if (not cf_f_act.is_sequence()) { cf_f_act.report_error( YAMLConfigError::ListExpected, "Fixed activity component must be provide as a list of" " triplet {label amount, unit}" ); } for (auto ind: RangeIterator(cf_f_act.size())) { auto cf = cf_f_act.get_section(ind); const auto label = cf.get_required_attribute< std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL); double value {0.0}; // value can be provided as value or log_value if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) { value = cf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT); } else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) { const auto log_val = cf.get_attribute< double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG); value = std::pow(10.0, log_val); } else { cf.report_error( YAMLConfigError::MissingRequiredAttribute, "Either amount or log amount must be provided." ); } reactant_box.add_fixed_activity_component(label, value); } } static void configure_constraints_fixed_fugacity( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_fug ) { if (not cf_f_fug.is_sequence()) { cf_f_fug.report_error( YAMLConfigError::ListExpected, "Fixed fugacity gas must be provide as a list of" " triplet {label amount, unit}" ); } for (auto ind: RangeIterator(cf_f_fug.size())) { auto cf = cf_f_fug.get_section(ind); // both component and gas must be provided const auto label_c = cf.get_required_attribute< std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT); const auto label_g = cf.get_required_attribute< std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL_GAS); double value {0.0}; // value can be provided as value or log_value if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) { value = cf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT); } else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) { const auto log_val = cf.get_attribute< double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG); value = std::pow(10.0, log_val); } else { cf.report_error( YAMLConfigError::MissingRequiredAttribute, "Either amount or log amount must be provided." ); } reactant_box.add_fixed_fugacity_gas(label_g, label_c, value); } } static void configure_constraints_fixed_molality( ReactantBox& reactant_box, YAMLConfigHandle&& cf_f_mol ) { if (not cf_f_mol.is_sequence()) { cf_f_mol.report_error( YAMLConfigError::ListExpected, "Fixed molality component must be provide as a list of" " triplet {label amount, unit}" ); } for (auto ind: RangeIterator(cf_f_mol.size())) { auto cf = cf_f_mol.get_section(ind); const auto label = cf.get_required_attribute< std::string>(SPC_CF_S_CONSTRAINTS_A_LABEL); double value {0.0}; // value can be provided as value or log_value if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) { value = cf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT); } else if (cf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) { const auto log_val = cf.get_attribute< double>(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG); value = std::pow(10.0, log_val); } else { cf.report_error( YAMLConfigError::MissingRequiredAttribute, "Either amount or log amount must be provided." ); } reactant_box.add_fixed_molality_component(label, value); } } } // end namespace internal // Constraints // =========== // new interface // ------------- void configure_specmicp_constraints( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints ) { // equation for water // ================== // enable conservation of water if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER)) { auto tmp_water = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER); if (tmp_water) { // only one constraints if (conf_constraints.get_optional_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM, false)) { conf_constraints.report_error(YAMLConfigError::InvalidArgument, "Attributes " SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "and " SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM " cannot be set to true at the same time"); } constraints.enable_conservation_water(); } else constraints.disable_conservation_water(); } // or enable saturated system if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) { if (conf_constraints.get_attribute(SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM)) { constraints.set_saturated_system(); } } conf_constraints.set_if_attribute_exists(constraints.inert_volume_fraction, SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION ); // Aqueous components // ================== // charge keeper // ------------- if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER)) { const auto label = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER); const index_t id = raw_db->get_id_component(label); if (id == no_species) { conf_constraints.report_error( YAMLConfigError::InvalidArgument, "Species " + label + " does not exist in the database." " It can't be the charge keeper." ); } constraints.set_charge_keeper(id); } // Fixed activity // --------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY)) { configure_specmicp_constraints_fixed_activity( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY) ); } // Fixed fugacity // -------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY)) { configure_specmicp_constraints_fixed_fugacity( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY) ); } // Fixed activity // --------------- if (conf_constraints.has_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY)) { configure_specmicp_constraints_fixed_molality( constraints, raw_db, conf_constraints.get_section(SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY) ); } // Surface model // ============= if (conf_constraints.has_attribute(SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL)) { const bool is_enabled = conf_constraints.get_attribute( SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL ); if (is_enabled) { const auto value = conf_constraints.get_required_attribute( SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION); constraints.enable_surface_model(value); } else { constraints.disable_surface_model(); } } } // helpers functions index_t get_id_component( const database::DataContainer* const raw_db, YAMLConfigHandle& conf, const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL ) { const auto label = conf.get_required_attribute(attr_key); const index_t id = raw_db->get_id_component(label); if (id == no_species) { conf.report_error( YAMLConfigError::InvalidArgument, "'" + label + "' is not a valid component." ); } return id; } index_t get_id_gas( const database::DataContainer* const raw_db, YAMLConfigHandle& conf, const std::string& attr_key=SPC_CF_S_CONSTRAINTS_A_LABEL_GAS ) { const auto label = conf.get_required_attribute(attr_key); const index_t id = raw_db->get_id_gas(label); if (id == no_species) { conf.report_error( YAMLConfigError::InvalidArgument, "'" + label + "' is not a valid gas." ); } return id; } scalar_t get_value(YAMLConfigHandle& conf) { scalar_t value = 0; if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)) { value = std::log10(conf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNT)); } else if (conf.has_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG)) { value = conf.get_attribute(SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG); } else { conf.report_error( YAMLConfigError::InvalidArgument, "Expected one of '" SPC_CF_S_CONSTRAINTS_A_AMOUNT "' or '" SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "'." ); } return value; } void configure_specmicp_constraints_fixed_activity( AdimensionalSystemConstraints& constraints, const database::DataContainer * const raw_db, YAMLConfigHandle&& conf_constraints_fixed_activity ) { for (uindex_t ind=0; ind(SPC_CF_S_INITIALIZATION_A_SOLUTION)); + } + if (conf_init.has_section(SPC_CF_S_INITIALIZATION_A_AQUEOUS)) + { + auto scf = conf_init.get_section(SPC_CF_S_INITIALIZATION_A_AQUEOUS); + if (not scf.is_map()) + { + scf.report_error(YAMLConfigError::MapExpected, + "A map {component, value} is expected."); + } + if (scf.has_attribute("all")) + { + solver.set_init_molality(scf.get_attribute("all")); + } + std::unordered_map tmap; + tmap.reserve(scf.size()); + for (auto it=scf.map_begin(); it!=scf.map_end(); ++it) + { + auto name = (*it).first; + if (name != "all") + { + tmap.insert({name, scf.get_attribute(name)}); + } + } + if (not tmap.empty()) + { + solver.set_init_molality(tmap); + } + } + if (conf_init.has_section(SPC_CF_S_INITIALIZATION_A_MINERALS)) + { + auto scf = conf_init.get_section(SPC_CF_S_INITIALIZATION_A_MINERALS); + if (not scf.is_map()) + { + scf.report_error(YAMLConfigError::MapExpected, + "A map {component, value} is expected."); + } + std::unordered_map tmap; + tmap.reserve(scf.size()); + for (auto it=scf.map_begin(); it!=scf.map_end(); ++it) + { + auto name = (*it).first; + tmap.insert({name, scf.get_attribute(name)}); + } + if (not tmap.empty()) + { + solver.set_init_volfrac_mineral(tmap); + } + } +} + } //end namespace io } //end namespace specmicp diff --git a/src/specmicp/io/configuration.hpp b/src/specmicp/io/configuration.hpp index c794177..955b3f2 100644 --- a/src/specmicp/io/configuration.hpp +++ b/src/specmicp/io/configuration.hpp @@ -1,102 +1,109 @@ /* ============================================================================= 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. * ============================================================================= */ #ifndef SPECMICP_SPECMICP_IO_HPP #define SPECMICP_SPECMICP_IO_HPP #include "specmicp_common/types.hpp" #include //! \file specmicp/io/configuration.hpp //! \brief configuration for specmicp adim system namespace specmicp { struct AdimensionalSystemSolverOptions; struct Formulation; struct AdimensionalSystemConstraints; class ReactantBox; +class SmartAdimSolver; namespace database { struct DataContainer; } using RawDatabasePtr = std::shared_ptr; namespace units { struct UnitsSet; } //end namespace units namespace io { class YAMLConfigHandle; /*! \brief Configure specmicp solver options Conserve previous values if not given. Configuration must be of the form : \code opts1: val1 opts2: val2 .... \endcode \warning units must be set elsewhere ! */ void SPECMICP_DLL_PUBLIC configure_specmicp_options( AdimensionalSystemSolverOptions& options, const units::UnitsSet& the_units, YAMLConfigHandle&& configuration ); //! \brief Configure a reactant box from the configuration ReactantBox SPECMICP_DLL_PUBLIC configure_specmicp_reactant_box( RawDatabasePtr raw_db, const units::UnitsSet& the_units, YAMLConfigHandle&& configuration ); //! \brief Configure specmicp constraints void SPECMICP_DLL_PUBLIC configure_specmicp_constraints( AdimensionalSystemConstraints& constraints, const database::DataContainer* const raw_db, YAMLConfigHandle&& conf_constraints ); +//! \brief Configure AdimSmartSolver initialization +void SPECMICP_DLL_PUBLIC configure_smart_solver_initialization( + SmartAdimSolver& solver, + YAMLConfigHandle&& conf_init + ); + } //end namespace io } //end namespace specmicp #endif // SPECMICP_SPECMICP_IO_HPP diff --git a/src/specmicp_common/io/config_yaml_sections.h b/src/specmicp_common/io/config_yaml_sections.h index 61917cc..1e320b2 100644 --- a/src/specmicp_common/io/config_yaml_sections.h +++ b/src/specmicp_common/io/config_yaml_sections.h @@ -1,325 +1,332 @@ /* ============================================================================= 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. * ============================================================================= */ #ifndef SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H #define SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H //! \file io/config_yaml_sections.h //! \brief Define sections in the YAML config files #define SPC_CF_NAME_PATH_ROOT "root" // common options // -------------- #define SPC_CF_A_MAX_ITER "maximum_iterations" #define SPC_CF_A_RES_TOL "residual_tolerance" #define SPC_CF_A_ABS_TOL "absolute_tolerance" #define SPC_CF_A_STEP_TOL "step_tolerance" #define SPC_CF_A_MAX_STEP_LENGTH "maximum_step_length" #define SPC_CF_A_MAX_STEP_MAX_ITER "maximum_step_maximum_iterations" #define SPC_CF_A_FILEPATH "file" #define SPC_CF_A_COMPONENT "component" #define SPC_CF_A_NODES "nodes" #define SPC_CF_A_TYPE "type" #define SPC_CF_V_DEFAULT "default" #define SPC_CF_V_PLUGIN "plugin" #define SPC_CF_A_PLUGIN_FILE "plugin_file" // filepath to a plugin #define SPC_CF_A_PLUGIN_NAME "plugin_name" // name of the plugin to use // user variables // -------------- #define SPC_CF_S_USRVARS "vars" #define SPC_CF_S_USRVARS_SS_INT "int" #define SPC_CF_S_USRVARS_SS_FLOAT "float" // plugin manager // -------------- #define SPC_CF_S_PLUGINS "plugins" #define SPC_CF_S_PLUGINS_A_SEARCHDIRS "dirs" #define SPC_CF_S_PLUGINS_A_TOLOAD "to_load" // Logs // ---- // type of output #define SPC_CF_S_LOGS_V_COUT "cout" #define SPC_CF_S_LOGS_V_CERR "cerr" #define SPC_CF_S_LOGS_V_FILE "file" // the different levels #define SPC_CF_S_LOGS_V_CRITICAL "critical" #define SPC_CF_S_LOGS_V_ERROR "error" #define SPC_CF_S_LOGS_V_WARNING "warning" #define SPC_CF_S_LOGS_V_INFO "info" #define SPC_CF_S_LOGS_V_DEBUG "debug" // runtime logs #define SPC_CF_S_LOGS "logs" #define SPC_CF_S_LOGS_A_LEVEL "level" #define SPC_CF_S_LOGS_A_OUTPUT "output" #define SPC_CF_S_LOGS_A_FILEPATH SPC_CF_A_FILEPATH // configuration logs #define SPC_CF_S_CONF_LOGS "configuration_logs" #define SPC_CF_S_CONF_LOGS_A_OUTPUT "output" #define SPC_CF_S_CONF_LOGS_A_FILE SPC_CF_A_FILEPATH // units // ----- #define SPC_CF_S_UNITS "units" #define SPC_CF_S_UNITS_A_LENGTH "length" #define SPC_CF_S_UNITS_A_MASS "mass" #define SPC_CF_S_UNITS_A_QUANTITY "quantity" // mesh // ---- #define SPC_CF_S_MESH "mesh" #define SPC_CF_S_MESH_A_TYPE "type" #define SPC_CF_S_MESH_SS_UNIFMESH "uniform_mesh" #define SPC_CF_S_MESH_SS_UNIFMESH_A_DX "dx" #define SPC_CF_S_MESH_SS_UNIFMESH_A_NBNODES "nb_nodes" #define SPC_CF_S_MESH_SS_UNIFMESH_A_SECTION "section" #define SPC_CF_S_MESH_SS_RAMPMESH "ramp_mesh" #define SPC_CF_S_MESH_SS_RAMPMESH_A_MIN_DX "min_dx" #define SPC_CF_S_MESH_SS_RAMPMESH_A_MAX_DX "max_dx" #define SPC_CF_S_MESH_SS_RAMPMESH_A_RAMP_LENGTH "length_ramp" #define SPC_CF_S_MESH_SS_RAMPMESH_A_PLATEAU_LENGTH "length_plateau" #define SPC_CF_S_MESH_SS_RAMPMESH_A_SECTION "section" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH "uniform_axisymmetric" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_RADIUS "radius" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_NBNODES "nb_nodes" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_HEIGHT "height" #define SPC_CF_S_MESH_SS_UNIFAXISYMESH_A_DX "dx" // database // -------- #define SPC_CF_S_DATABASE "database" #define SPC_CF_S_DATABASE_A_PATH "path" #define SPC_CF_S_DATABASE_A_SWAP "swap_components" #define SPC_CF_S_DATABASE_A_SWAP_K_IN "in" #define SPC_CF_S_DATABASE_A_SWAP_K_OUT "out" #define SPC_CF_S_DATABASE_A_REMOVE_GAS "remove_gas" #define SPC_CF_S_DATABASE_A_REMOVE_SOLIDS "remove_solids" #define SPC_CF_S_DATABASE_A_REMOVE_SORBED "remove_sorbeds" #define SPC_CF_S_DATABASE_A_EXTRA_GAS "extra_gas" #define SPC_CF_S_DATABASE_A_EXTRA_SOLIDS "extra_solids" #define SPC_CF_S_DATABASE_A_EXTRA_SORBED "extra_sorbeds" #define SPC_CF_S_DATABASE_A_REMOVE_HALF_CELLS "remove_half_cells" #define SPC_CF_S_DATABASE_A_LIST_SOLIDS_TOKEEP "list_solids_to_keep" // SpecMiCP options // ---------------- #define SPC_CF_S_SPECMICP "specmicp_options" #define SPC_CF_S_SPECMICP_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_SPECMICP_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_SPECMICP_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_SPECMICP_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH #define SPC_CF_S_SPECMICP_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER #define SPC_CF_S_SPECMICP_A_SCALING "enable_scaling" #define SPC_CF_S_SPECMICP_A_NONMONOTONE "enable_nonmonotone_linesearch" #define SPC_CF_S_SPECMICP_A_DESCENT_DIRECTION "factor_descent_direction" #define SPC_CF_S_SPECMICP_A_COND_CHECK "threshold_condition_check" #define SPC_CF_S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH "threshold_cycling_linesearch" #define SPC_CF_S_SPECMICP_A_NONIDEAL_TOL "non_ideality_tolerance" #define SPC_CF_S_SPECMICP_A_NONIDEAL_MAX_ITER "non_ideality_" SPC_CF_A_MAX_ITER #define SPC_CF_S_SPECMICP_A_CUTOFF_TOT_CONC "cutoff_total_concentration" #define SPC_CF_S_SPECMICP_A_RESTART_CONCENTRATION "restart_concentration" #define SPC_CF_S_SPECMICP_A_RESTART_WATER_VOL_FRAC "restart_water_volume_fraction" #define SPC_CF_S_SPECMICP_A_UNDER_RELAXATION "under_relaxation" // SpecMiCP formulation // -------------------- #define SPC_CF_S_FORMULATION "formulation" #define SPC_CF_S_FORMULATION_A_SOLUTION "solution" #define SPC_CF_S_FORMULATION_A_AQUEOUS "aqueous" #define SPC_CF_S_FORMULATION_A_MINERALS "solid_phases" #define SPC_CF_S_FORMULATION_A_AMOUNT "amount" #define SPC_CF_S_FORMULATION_A_UNIT "unit" #define SPC_CF_S_FORMULATION_A_LABEL "label" +#define SPC_CF_S_INITIALIZATION "initialization" +#define SPC_CF_S_INITIALIZATION_A_SOLUTION "solution" +#define SPC_CF_S_INITIALIZATION_A_MINERALS "solid_phases" +#define SPC_CF_S_INITIALIZATION_A_AQUEOUS "aqueous" + + + // SpecMiCP constraints // -------------------- #define SPC_CF_S_CONSTRAINTS "constraints" #define SPC_CF_S_CONSTRAINTS_A_CHARGEKEEPER "charge_keeper" #define SPC_CF_S_CONSTRAINTS_A_FIXEDACTIVITY "fixed_activity" #define SPC_CF_S_CONSTRAINTS_A_FIXEDMOLALITY "fixed_molality" #define SPC_CF_S_CONSTRAINTS_A_FIXEDFUGACITY "fixed_fugacity" #define SPC_CF_S_CONSTRAINTS_A_SATURATED_SYSTEM "saturated_system" #define SPC_CF_S_CONSTRAINTS_A_CONSERVATION_WATER "conservation_water" #define SPC_CF_S_CONSTRAINTS_A_SURFACE_MODEL "surface_model" #define SPC_CF_S_CONSTRAINTS_A_SURFACE_CONCENTRATION "surface_site_concentration" #define SPC_CF_S_CONSTRAINTS_A_INERT_VOLUME_FRACTION "inert_volume_fraction" #define SPC_CF_S_CONSTRAINTS_A_LABEL "label" #define SPC_CF_S_CONSTRAINTS_A_AMOUNT "amount" #define SPC_CF_S_CONSTRAINTS_A_AMOUNTLOG "amount_log" #define SPC_CF_S_CONSTRAINTS_A_LABEL_COMPONENT "label_component" #define SPC_CF_S_CONSTRAINTS_A_LABEL_GAS "label_gas" // Variables // --------- #define SPC_CF_S_REACTMICPVARIABLES "variables" #define SPC_CF_S_REACTMICPVARIABLES_A_TYPE SPC_CF_A_TYPE // Boundary conditions // ------------------- #define SPC_CF_S_BOUNDARYCONDITIONS "boundary_conditions" // transport options // ----------------- #define SPC_CF_S_DFPM "transport_options" #define SPC_CF_S_DFPM_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_DFPM_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_DFPM_A_ABS_TOL SPC_CF_A_ABS_TOL #define SPC_CF_S_DFPM_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_DFPM_A_TRSHOLD_STATIONARY "threshold_stationary" #define SPC_CF_S_DFPM_A_MAX_STEP_LENGTH SPC_CF_A_MAX_STEP_LENGTH #define SPC_CF_S_DFPM_A_MAX_STEP_MAX_ITER SPC_CF_A_MAX_STEP_MAX_ITER #define SPC_CF_S_DFPM_A_SPARSE_SOLVER "sparse_solver" #define SPC_CF_S_DFPM_A_LINESEARCH "linesearch" #define SPC_CF_S_DFPM_A_QUASI_NEWTON "quasi_newton" // ReactMiCP coupling solver // ------------------------- #define SPC_CF_S_REACTMICP "reactmicp_options" #define SPC_CF_S_REACTMICP_A_RES_TOL SPC_CF_A_RES_TOL #define SPC_CF_S_REACTMICP_A_ABS_TOL SPC_CF_A_ABS_TOL #define SPC_CF_S_REACTMICP_A_STEP_TOL SPC_CF_A_STEP_TOL #define SPC_CF_S_REACTMICP_A_GOOD_ENOUGH_TOL "good_enough_tolerance" #define SPC_CF_S_REACTMICP_A_MAX_ITER SPC_CF_A_MAX_ITER #define SPC_CF_S_REACTMICP_A_IMPL_UPSCALING "implicit_upscaling" // ReactMiCP output // ---------------- #define SPC_CF_S_REACTOUTPUT "reactmicp_output" #define SPC_CF_S_REACTOUTPUT_A_FILEPATH SPC_CF_A_FILEPATH #define SPC_CF_S_REACTOUTPUT_A_TYPE "type" #define SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5 "hdf5" #define SPC_CF_S_REACTOUTPUT_A_TYPE_V_CSV "csv" #define SPC_CF_S_REACTOUTPUT_A_POROSITY "porosity" #define SPC_CF_S_REACTOUTPUT_A_PH "ph" #define SPC_CF_S_REACTOUTPUT_A_DIFFUSIVITY "diffusivity" #define SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL "volume_fraction_solid" #define SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC "total_aqueous_concentration" #define SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC "total_solid_concentration" #define SPC_CF_S_REACTOUTPUT_A_TOT_CONC "total_concentration" #define SPC_CF_S_REACTOUTPUT_A_SATINDEX "saturation_index" #define SPC_CF_S_REACTOUTPUT_A_DATABASE "database" // ReactMiCP simulation stuff // -------------------------- #define SPC_CF_S_SIMULINFO "simulation" #define SPC_CF_S_SIMULINFO_A_NAME "name" #define SPC_CF_S_SIMULINFO_A_OUTPUTPREFIX "output_prefix" #define SPC_CF_S_SIMULINFO_A_PRINTITER "print_iter_info" #define SPC_CF_S_SIMULINFO_A_OUTPUTSTEP "output_step" // ReactMiCP adaptive timestep // --------------------------- #define SPC_CF_S_TIMESTEPPER "timestepper" #define SPC_CF_S_TIMESTEP_A_LOWER_BOUND "minimum_dt" #define SPC_CF_S_TIMESTEP_A_UPPER_BOUND "maximum_dt" #define SPC_CF_S_TIMESTEP_A_RESTART_DT "restart_dt" #define SPC_CF_S_TIMESTEP_A_LOWER_TARGET "lower_iterations_target" #define SPC_CF_S_TIMESTEP_A_UPPER_TARGET "upper_iterations_target" #define SPC_CF_S_TIMESTEP_A_AVERAGE_PARAMETER "average_parameter" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_FAILURE "decrease_factor_if_failure" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_MINIMUM "increase_factor_if_minumum" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_DECREASE "decrease_factor" #define SPC_CF_S_TIMESTEP_A_FACTOR_IF_INCREASE "increase_factor" // ReactMiCP Staggers // ------------------ #define SPC_CF_S_STAGGERS "staggers" #define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER "upscaling" #define SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE SPC_CF_A_TYPE #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER "chemistry" #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE SPC_CF_A_TYPE #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM "equilibrium" #define SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS "equilibrium_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER "transport" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION "merge_saturation" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS "merge_aqueous" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_CUTOFFRESIDUAL "cutoff_R_0_squared" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS "default_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS "saturation_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS "aqueous_options" #define SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS "gas_options" // Run ! // ----- #define SPC_CF_S_RUN "run" #define SPC_CF_S_RUN_A_RUNUTNIL "run_until" #endif // SPECMICP_UTILS_IO_CONFIGYAMLSECTIONS_H diff --git a/src/specmicp_common/io/safe_config.cpp b/src/specmicp_common/io/safe_config.cpp index a599919..8e8a37b 100644 --- a/src/specmicp_common/io/safe_config.cpp +++ b/src/specmicp_common/io/safe_config.cpp @@ -1,810 +1,819 @@ /* ============================================================================= 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 "safe_config.hpp" #include "specmicp_common/string_algorithms.hpp" #include #include #include #include #include "specmicp_common/log.hpp" #include "specmicp_common/compat.hpp" #define SECTION_USER_VARIABLES "vars" #define SECTION_INT_USER_VARIABLES "int" #define SECTION_FLOAT_USER_VARIABLES "float" #define NAME_PATH_ROOT "root" namespace specmicp { namespace io { // ================ // ConfigFileHandle // ================ struct YAMLConfigFileHandle { std::string m_filename; std::unordered_map m_float_vars; std::unordered_map m_int_vars; YAMLConfigFileHandle(const std::string& filename): m_filename(filename) {} }; // ============ // ConfigHandle // ============ struct YAMLConfigHandle::YAMLConfigHandleImpl { YAML::Node m_node; std::shared_ptr m_file; std::string m_path; std::string m_section_name; std::vector m_asked_attrs; YAMLConfigHandleImpl( const YAML::Node& node, std::shared_ptr const file_handle, const std::string& section, const std::string& path): m_node(node), m_file(file_handle), m_path(path), m_section_name(section) { if (node.IsMap()) { m_asked_attrs.reserve(node.size()); } } ~YAMLConfigHandleImpl(); scalar_t parse_scalar_expr(const std::string& expr) { return utils::parse_expression( expr, m_file->m_float_vars ); } index_t parse_int_expr(const std::string& expr) { return utils::parse_expression( expr, m_file->m_int_vars ); } bool parse_bool_expr(const std::string& expr) { return utils::string_to_bool(expr); } void record_attr(const std::string& attr) { m_asked_attrs.push_back(attr); } void assert_open() { if (m_file == nullptr) { throw std::runtime_error( "Information from file '" + m_file->m_filename + "' are lost..."); } } template std::vector list_to_vector( const std::string& list, YAMLConfigHandle& parent ); }; YAMLConfigHandle::YAMLConfigHandle(const YAML::Node& node, const std::shared_ptr file_handle, const std::string& section, const std::string& path): m_impl(utils::make_pimpl( node, file_handle, section, path)) { } YAMLConfigHandle& YAMLConfigHandle::operator= (const YAMLConfigHandle& other) { m_impl = other.m_impl; m_impl->assert_open(); return *this; } YAMLConfigHandle::YAMLConfigHandle(const YAMLConfigHandle& other): m_impl(other.m_impl) { m_impl->assert_open(); } YAMLConfigHandle::YAMLConfigHandle(YAMLConfigHandle&& other): m_impl(std::move(other.m_impl)) { m_impl->assert_open(); } YAMLConfigHandle::~YAMLConfigHandle() = default; uindex_t YAMLConfigHandle::size() { return static_cast(m_impl->m_node.size()); } void YAMLConfigHandle::set_file_handle( std::shared_ptr file_handle ) { m_impl->m_file = file_handle; m_impl->assert_open(); } bool YAMLConfigHandle::is_map() { return m_impl->m_node.IsMap(); } bool YAMLConfigHandle::is_map(const std::string& node) { bool val = false; if (m_impl->m_node[node]) { val = m_impl->m_node[node].IsMap(); } return val; } bool YAMLConfigHandle::is_sequence() { return m_impl->m_node.IsSequence(); } bool YAMLConfigHandle::is_sequence(const std::string& node) { bool val = false; if (m_impl->m_node[node]) { val = m_impl->m_node[node].IsSequence(); } return val; } bool YAMLConfigHandle::has_node(const std::string& node) { return (m_impl->m_node[node]); } bool YAMLConfigHandle::has_node(uindex_t index) { return (m_impl->m_node[index]); } bool YAMLConfigHandle::has_section(const std::string& section) { bool val = false; if (m_impl->m_node[section]) { YAML::Node attr = m_impl->m_node[section]; val = (attr.IsMap() or attr.IsSequence()); } return val; } bool YAMLConfigHandle::has_section(uindex_t index) { bool val = false; if (m_impl->m_node[index]) { YAML::Node attr = m_impl->m_node[index]; val = (attr.IsMap() or attr.IsSequence()); } return val; } YAMLConfigHandle YAMLConfigHandle::get_section(const std::string& section) { if (not has_section(section)) { report_error(YAMLConfigError::MissingRequiredSection, section); } m_impl->record_attr(section); return YAMLConfigHandle(m_impl->m_node[section], m_impl->m_file, section, m_impl->m_path+"->"+section); } YAMLConfigHandle YAMLConfigHandle::get_section(uindex_t index) { auto section_name = "["+std::to_string(index)+"]"; if (not has_section(index)) { report_error(YAMLConfigError::MissingRequiredSection, section_name ); } return YAMLConfigHandle(m_impl->m_node[index], m_impl->m_file, section_name, m_impl->m_path+section_name); } bool YAMLConfigHandle::has_attribute(const std::string& attribute) { bool val = false; if (m_impl->m_node[attribute]) { YAML::Node attr = m_impl->m_node[attribute]; val = not (attr.IsMap() or attr.IsSequence()); } return val; } #define report_conversion_error(attribute, type) \ report_error( \ YAMLConfigError::ConversionError, \ "Conversion of attribute '" \ + attribute + "' to " type + " : " + e.what() \ ); template <> scalar_t YAMLConfigHandle::get_attribute(const std::string& attribute) { m_impl->record_attr(attribute); YAML::Node node = m_impl->m_node[attribute]; scalar_t val; m_impl->assert_open(); try { val = m_impl->parse_scalar_expr(node.as()); } catch (const std::invalid_argument& e) { report_conversion_error(attribute, "scalar"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of attribute '" + attribute + "' to scalar : " + e.what() ); } return val; } template <> index_t YAMLConfigHandle::get_attribute(const std::string& attribute) { m_impl->record_attr(attribute); YAML::Node node = m_impl->m_node[attribute]; index_t val; m_impl->assert_open(); try { val = m_impl->parse_int_expr(node.as()); } catch (const std::invalid_argument& e) { report_conversion_error(attribute, "integer"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of attribute '" + attribute + "' to integer : " + e.what() ); } return val; } template <> uindex_t YAMLConfigHandle::get_attribute(const std::string& attribute) { m_impl->record_attr(attribute); YAML::Node node = m_impl->m_node[attribute]; index_t val; m_impl->assert_open(); try { val = m_impl->parse_int_expr(node.as()); } catch (const std::invalid_argument& e) { report_conversion_error(attribute, "integer"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of attribute '" + attribute + "' to integer : " + e.what() ); } if (val < 0) { report_error(YAMLConfigError::InvalidArgument, "Expected positive value for attribute " + attribute + ", got '" + std::to_string(val) + "'."); } return val; } template <> bool YAMLConfigHandle::get_attribute(const std::string& attribute) { m_impl->record_attr(attribute); YAML::Node node = m_impl->m_node[attribute]; bool val; try { val = m_impl->parse_bool_expr(node.as()); } catch (const std::invalid_argument& e) { report_conversion_error(attribute, "boolean"); } return val; } template <> std::string YAMLConfigHandle::get_attribute(const std::string& attribute) { m_impl->record_attr(attribute); return m_impl->m_node[attribute].as(); } #undef report_converssion_error // list #define report_list_conversion_error(index, type) \ report_error( \ YAMLConfigError::ConversionError, \ "Conversion of value at index '" \ + std::to_string(index) + "' to " \ type + " : " + e.what() \ ); template <> scalar_t YAMLConfigHandle::get_value(uindex_t ind) { YAML::Node node = m_impl->m_node[ind]; scalar_t val; m_impl->assert_open(); try { val = m_impl->parse_scalar_expr(node.as()); } catch (const std::invalid_argument& e) { report_list_conversion_error(ind, "scalar"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of value at index '" + std::to_string(ind) + "' to scalar : " + e.what() ); } return val; } template <> index_t YAMLConfigHandle::get_value(uindex_t ind) { YAML::Node node = m_impl->m_node[ind]; index_t val; m_impl->assert_open(); try { val = m_impl->parse_int_expr(node.as()); } catch (const std::invalid_argument& e) { report_list_conversion_error(ind, "integer"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of value at index '" + std::to_string(ind) + "' to integer : " + e.what() ); } return val; } template <> uindex_t YAMLConfigHandle::get_value(uindex_t ind) { YAML::Node node = m_impl->m_node[ind]; index_t val; m_impl->assert_open(); try { val = m_impl->parse_int_expr(node.as()); } catch (const std::invalid_argument& e) { report_list_conversion_error(ind, "integer"); } catch (const std::out_of_range& e) { report_error(YAMLConfigError::UnknownVariable, "Conversion of value at index '" + std::to_string(ind) + "' to integer : " + e.what() ); } if (val < 0) { report_error(YAMLConfigError::InvalidArgument, "Expected positive value at index " + std::to_string(ind) + ", got '" + std::to_string(val) + "'."); } return val; } template <> bool YAMLConfigHandle::get_value(uindex_t ind) { YAML::Node node = m_impl->m_node[ind]; bool val; try { val = m_impl->parse_bool_expr(node.as()); } catch (const std::invalid_argument& e) { report_list_conversion_error(ind, "boolean"); } return val; } template <> std::string YAMLConfigHandle::get_value(uindex_t ind) { return m_impl->m_node[ind].as(); } void YAMLConfigHandle::report_error( YAMLConfigError error_type, const std::string& error_msg ) { std::string header; m_impl->assert_open(); if (m_impl->m_file->m_filename != "") { header += "\n\t - in file " + m_impl->m_file->m_filename; } if (m_impl->m_path != "") { header += "\n\t - in section " + m_impl->m_path + "\n"; } switch (error_type) { case YAMLConfigError::ConversionError: throw std::invalid_argument("Conversion error : " + header + error_msg); break; case YAMLConfigError::UnknownVariable: throw std::out_of_range("Unknown variable :" + header + error_msg); break; case YAMLConfigError::MissingRequiredAttribute: throw std::out_of_range("Missing required attribute : " + error_msg + header); break; case YAMLConfigError::MissingRequiredSection: throw std::out_of_range("Missing required section : " + error_msg + header); break; default: throw std::runtime_error(error_msg + header); break; } } YAMLConfigHandle::YAMLConfigHandleImpl::~YAMLConfigHandleImpl() { if (m_node.IsMap()) { assert_open(); // report untouched attribute std::vector unread; for (auto it=m_node.begin(); it!=m_node.end(); ++it) { auto is_read = std::find( m_asked_attrs.cbegin(), m_asked_attrs.cend(), it->first.as() ); if (is_read == m_asked_attrs.cend()) { WARNING << "Unread key : " + it->first.as() + " in file '" + m_file->m_filename + "'" + " in section '" + m_path + "'."; } } } } // =========== // MapIterator // =========== +YAMLConfigHandle::MapIterator::MapIterator(MapIterator&& other): + m_handle(other.m_handle), + m_true_it(other.m_true_it.release()) +{ + +} + +YAMLConfigHandle::MapIterator::~MapIterator() = default; + std::pair YAMLConfigHandle::MapIterator::operator* () { m_handle->record_attr((*m_true_it)->first.as()); return { (*m_true_it)->first.as(), (*m_true_it)->second.as() }; } YAMLConfigHandle::MapIterator& YAMLConfigHandle::MapIterator::operator++ () { m_true_it->operator++ (); return *this; } bool YAMLConfigHandle::MapIterator::operator==(const MapIterator& other) { return (*m_true_it == *other.m_true_it); } bool YAMLConfigHandle::MapIterator::operator!=(const MapIterator& other) { return (*m_true_it != *other.m_true_it); } YAMLConfigHandle::MapIterator::MapIterator( YAMLConfigHandleImpl* handle, std::unique_ptr it ): m_handle(handle), m_true_it(it.release()) {} YAMLConfigHandle::MapIterator YAMLConfigHandle::map_begin() { return MapIterator( m_impl.get(), make_unique(m_impl->m_node.begin()) ); } YAMLConfigHandle::MapIterator YAMLConfigHandle::map_end() { return MapIterator( m_impl.get(), make_unique(m_impl->m_node.end()) ); } // List to vector // ------------- // implementation of list_to_vector template std::vector YAMLConfigHandle::YAMLConfigHandleImpl::list_to_vector( const std::string& list, YAMLConfigHandle& parent ) { // check that it is a list as expected if (not m_node[list]) { parent.report_error(YAMLConfigError::MissingRequiredAttribute, list); } auto node = m_node[list]; if (not node.IsSequence()) { parent.report_error(YAMLConfigError::ListExpected, "A list was expected for attribute : " + list); } record_attr(list); // parse the list std::vector vec; vec.reserve(parent.size()); auto list_section = parent.get_section(list); uindex_t size = list_section.size(); for (uindex_t ind=0; ind(ind)); } return vec; } // two specialization for indices => also read range of indices template <> std::vector YAMLConfigHandle::YAMLConfigHandleImpl::list_to_vector( const std::string& list, YAMLConfigHandle& parent ) { if (not m_node[list]) { parent.report_error(YAMLConfigError::MissingRequiredAttribute, list); } record_attr(list); // parse the list auto node = m_node[list]; std::vector vec; if (not node.IsSequence()) { try { // it ;ay be our special way to generate list "1:10" vec = utils::range_indices(node.as()); } catch (const std::invalid_argument& e) { parent.report_error(YAMLConfigError::InvalidArgument, "Expected valid list for attribute " + list + ". Error : " + e.what()); } } else { vec.reserve(parent.size()); auto list_section = parent.get_section(list); uindex_t size = list_section.size(); for (uindex_t ind=0; ind(ind)); } } return vec; } template <> std::vector YAMLConfigHandle::YAMLConfigHandleImpl::list_to_vector( const std::string& list, YAMLConfigHandle& parent ) { if (not m_node[list]) { parent.report_error(YAMLConfigError::MissingRequiredAttribute, list); } record_attr(list); // parse the list auto node = m_node[list]; std::vector vec; if (not node.IsSequence()) { try { // it may be the slice way to generate lists ":" vec = utils::range_indices(node.as()); } catch (const std::invalid_argument& e) { parent.report_error(YAMLConfigError::InvalidArgument, "Expected valid list for attribute " + list + ". Error : " + e.what()); } } else { vec.reserve(parent.size()); auto list_section = parent.get_section(list); uindex_t size = list_section.size(); for (uindex_t ind=0; ind(ind)); } } return vec; } template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list) { return m_impl->list_to_vector(list, *this); } template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list) { return m_impl->list_to_vector(list, *this); } template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list) { return m_impl->list_to_vector(list, *this); } template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list) { return m_impl->list_to_vector(list, *this); } template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list) { return m_impl->list_to_vector(list, *this); } // ========== // ConfigFile // ========== YAMLConfigFile YAMLConfigFile::load( const std::string& file ) { return YAMLConfigFile(YAML::LoadFile(file), file); } YAMLConfigFile YAMLConfigFile::load_from_string( const std::string& str, std::string name ) { return YAMLConfigFile(YAML::Load(str), name); } std::unique_ptr YAMLConfigFile::make( const std::string& file ) { // can't use make unique, YAMLConfigFile() is private.... return std::unique_ptr( new YAMLConfigFile(YAML::LoadFile(file), file) ); } std::unique_ptr YAMLConfigFile::make_from_string( const std::string& str, std::string name ) { return std::unique_ptr( new YAMLConfigFile(YAML::Load(str), name) ); } YAMLConfigFile::YAMLConfigFile( const YAML::Node& node, const std::string name ): YAMLConfigHandle(node, nullptr, NAME_PATH_ROOT, NAME_PATH_ROOT), m_handle(std::make_shared(name)) { set_file_handle(m_handle); // read variables if (has_section(SECTION_USER_VARIABLES)) { auto& int_vars = m_handle->m_int_vars; auto& float_vars = m_handle->m_float_vars; auto var_section = get_section(SECTION_USER_VARIABLES); if (var_section.has_section(SECTION_INT_USER_VARIABLES)) { auto int_section = var_section.get_section(SECTION_INT_USER_VARIABLES); int_vars.reserve( int_vars.size() + int_section.size()); float_vars.reserve(float_vars.size() + int_section.size()); for ( auto it=int_section.map_begin(); it!=int_section.map_end(); ++it ) { auto tp = *it; auto key = utils::strip(tp.first); index_t val = utils::parse_expression(tp.second, int_vars); int_vars.insert({key, val}); float_vars.insert({key, static_cast(val)}); } } if (var_section.has_section(SECTION_FLOAT_USER_VARIABLES)) { auto float_section = var_section.get_section(SECTION_FLOAT_USER_VARIABLES); float_vars.reserve(float_vars.size()+float_section.size()); for ( auto it=float_section.map_begin(); it!=float_section.map_end(); ++it ) { auto tp = *it; auto key = utils::strip(tp.first); scalar_t val = utils::parse_expression(tp.second); float_vars.insert({key, val}); } } } } YAMLConfigFile::~YAMLConfigFile() = default; } //end namespace io } //end namespace specmicp diff --git a/src/specmicp_common/io/safe_config.hpp b/src/specmicp_common/io/safe_config.hpp index 2ccc539..1ac825b 100644 --- a/src/specmicp_common/io/safe_config.hpp +++ b/src/specmicp_common/io/safe_config.hpp @@ -1,530 +1,534 @@ /* ============================================================================= 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. * ============================================================================= */ #ifndef SPECMICP_UTILS_SAFECONFIG_HPP #define SPECMICP_UTILS_SAFECONFIG_HPP //! \file safe_config.hpp //! \brief YAML configuration file reader //! //! This is a set of wrappers over the yaml-cpp API #include "specmicp_common/types.hpp" #include "specmicp_common/pimpl_ptr.hpp" #include #include #include namespace YAML { class Node; } //end namespace YAML namespace specmicp { namespace io { //! \brief Type of error that can be reported enum class YAMLConfigError { UnknownError, UnknownVariable, MissingRequiredAttribute, MissingRequiredSection, InvalidArgument, ConversionError, ListExpected, MapExpected }; // internal store required information //! \internal struct SPECMICP_DLL_LOCAL YAMLConfigFileHandle; //! \brief Access to the configuration class SPECMICP_DLL_PUBLIC YAMLConfigHandle { public: //! \brief Report an error void report_error( YAMLConfigError error_type, const std::string& error_msg ); //! \brief Return true if it is a sequence bool is_sequence(); //! \brief Return true if sub-node exists and is a sequence bool is_sequence(const std::string& node); //! \brief Return true if sub-node exists and is a map bool is_map(const std::string& node); //! \brief Return true if it is a map bool is_map(); //! \brief Return true if sub-node exist bool has_node(const std::string& node); //! \brief Return true if sub-node exist bool has_node(uindex_t value); //! \brief Return true if the node has the given attribute //! //! An attribute is a pair key/value, where the value is //! neither a map or a sequence bool has_attribute(const std::string& attribute); //! \brief Return true if the node has the given section //! //! A section is a pair key/value, where the value is a //! map or a sequence bool has_section(const std::string& section); //! \brief Return true if the node has the given section //! //! A section is a pair key/value, where the value is a //! map or a sequence bool has_section(uindex_t value); //! \brief Return the number of element in the section uindex_t size(); //! \brief Return a section from a map YAMLConfigHandle get_section(const std::string& section); //! \brief Return a section from a list YAMLConfigHandle get_section(uindex_t value); //! \brief Return an attribute template T get_attribute(const std::string& attribute); //! \brief Return a lower bounded scalar attribute template T get_attribute(const std::string& attribute, T lower_bound ); //! \brief Return a bounded scalar attribute template T get_attribute(const std::string& attribute, T lower_bound, T upper_bound ); //! \brief Return a value in a list template T get_value(uindex_t ind); //! \brief Return a required attribute template T get_required_attribute(const std::string& attribute); //! \brief Return a required attribute template T get_required_attribute( const std::string& attribute, T lower_bound ); //! \brief Return a required attribute template T get_required_attribute( const std::string& attribute, T lower_bound, T upper_bound ); //! \brief Return an optional attribute template T get_optional_attribute( const std::string& attribute, const T& default_value ); //! \brief Return an optional attribute template T get_optional_attribute( const std::string& attribute, const T& default_value, T lower_bound ); //! \brief Return an optional attribute template T get_optional_attribute( const std::string& attribute, const T& default_value, T lower_bound, T upper_bound ); //! \brief Set a variable if the attribute exists template void set_if_attribute_exists( T& var_to_set, const std::string& attribute ); template void set_if_attribute_exists( T& var_to_set, const std::string& attribute, T lower_bound ); template void set_if_attribute_exists( T& var_to_set, const std::string& attribute, T lower_bound, T upper_bound ); template void set_if_attribute_exists( int& var_to_set, const std::string& attribute ); template void set_if_attribute_exists( int& var_to_set, const std::string& attribute, T lower_bound ); //! \brief Return a list of values for a list attribute template std::vector list_to_vector(const std::string& list); ~YAMLConfigHandle(); YAMLConfigHandle& operator= (const YAMLConfigHandle& other); YAMLConfigHandle(const YAMLConfigHandle& other); YAMLConfigHandle(YAMLConfigHandle&& other); protected: // can only be created from a file or another handle YAMLConfigHandle( const YAML::Node& node, std::shared_ptr const file_handle, const std::string& section, const std::string& path); void set_file_handle(std::shared_ptr file_handle); private: struct SPECMICP_DLL_LOCAL YAMLConfigHandleImpl; utils::pimpl_ptr m_impl; public: //! Iterator over a map class MapIterator { friend YAMLConfigHandle; public: + MapIterator(MapIterator&& other); + ~MapIterator(); + std::pair operator* (); + MapIterator& operator++ (); bool operator==(const MapIterator& other); bool operator!=(const MapIterator& other); private: MapIterator( YAMLConfigHandleImpl* handle, std::unique_ptr it ); YAMLConfigHandleImpl* m_handle; std::unique_ptr m_true_it; }; MapIterator map_begin(); MapIterator map_end(); }; // ======================= // // YAMLConfigFile // // ======================= // //! \brief This class represent a file //! //! It must be alive while the config is accessed class SPECMICP_DLL_PUBLIC YAMLConfigFile: public YAMLConfigHandle { public: static YAMLConfigFile load( const std::string& file ); static YAMLConfigFile load_from_string( const std::string& str, std::string name="" ); static std::unique_ptr make( const std::string& file ); static std::unique_ptr make_from_string( const std::string& str, std::string name="" ); ~YAMLConfigFile(); private: YAMLConfigFile( const YAML::Node& node, const std::string name ); std::shared_ptr m_handle; }; // ======================= // // Template implementation // // ======================= // // explicit specialization template <> scalar_t YAMLConfigHandle::get_attribute(const std::string& attribute); template <> index_t YAMLConfigHandle::get_attribute(const std::string& attribute); template <> uindex_t YAMLConfigHandle::get_attribute(const std::string& attribute); template <> bool YAMLConfigHandle::get_attribute(const std::string& attribute); template <> std::string YAMLConfigHandle::get_attribute(const std::string& attribute); template T YAMLConfigHandle::get_attribute( const std::string& attribute, T lower_bound, T upper_bound) { T value = get_attribute(attribute); if (value < lower_bound or value > upper_bound) { report_error( YAMLConfigError::InvalidArgument, "Attribute '" + attribute + "' is required to" " be between " + std::to_string(lower_bound) + " and " + std::to_string(upper_bound) + " (Got : " + std::to_string(value) + ")." ); } return value; } template T YAMLConfigHandle::get_attribute( const std::string& attribute, T lower_bound ) { T value = get_attribute(attribute); if (value < lower_bound ) { report_error( YAMLConfigError::InvalidArgument, "Attribute '" + attribute + "' is required to" " be greater than " + std::to_string(lower_bound) + " (Got : " + std::to_string(value) + ")." ); } return value; } template <> scalar_t YAMLConfigHandle::get_value(uindex_t ind); template <> index_t YAMLConfigHandle::get_value(uindex_t ind); template <> uindex_t YAMLConfigHandle::get_value(uindex_t ind); template <> bool YAMLConfigHandle::get_value(uindex_t ind); template <> std::string YAMLConfigHandle::get_value(uindex_t ind); template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list); template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list); template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list); template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list); template <> std::vector YAMLConfigHandle::list_to_vector(const std::string& list); // Required attribute // ------------------ template T YAMLConfigHandle::get_required_attribute(const std::string& attribute) { if (not has_attribute(attribute)) { report_error(YAMLConfigError::MissingRequiredAttribute, attribute); } return get_attribute(attribute); } template T YAMLConfigHandle::get_required_attribute( const std::string& attribute, T lower_bound ) { if (not has_attribute(attribute)) { report_error(YAMLConfigError::MissingRequiredAttribute, attribute); } return get_attribute(attribute, lower_bound); } template T YAMLConfigHandle::get_required_attribute( const std::string& attribute, T lower_bound, T upper_bound ) { if (not has_attribute(attribute)) { report_error(YAMLConfigError::MissingRequiredAttribute, attribute); } return get_attribute(attribute, lower_bound, upper_bound); } // Optional attribute // ------------------ template T YAMLConfigHandle::get_optional_attribute( const std::string& attribute, const T& default_value) { if (not has_attribute(attribute)) { return default_value; } return get_attribute(attribute); } template T YAMLConfigHandle::get_optional_attribute( const std::string& attribute, const T& default_value, T lower_bound ) { if (not has_attribute(attribute)) { return default_value; } return get_attribute(attribute, lower_bound); } //! \brief Return an optional attribute template T YAMLConfigHandle::get_optional_attribute( const std::string& attribute, const T& default_value, T lower_bound, T upper_bound ) { if (not has_attribute(attribute)) { return default_value; } return get_attribute(attribute, lower_bound, upper_bound); } // Set if attribute exists // ----------------------- template void YAMLConfigHandle::set_if_attribute_exists( T& var_to_set, const std::string& attribute ) { if (has_attribute(attribute)) { var_to_set = get_attribute(attribute); } } template void YAMLConfigHandle::set_if_attribute_exists( int& var_to_set, const std::string& attribute ) { if (has_attribute(attribute)) { var_to_set = static_cast(get_attribute(attribute)); } } template void YAMLConfigHandle::set_if_attribute_exists( T& var_to_set, const std::string& attribute, T lower_bound, T upper_bound ) { if (has_attribute(attribute)) { var_to_set = get_attribute(attribute, lower_bound, upper_bound); } } template void YAMLConfigHandle::set_if_attribute_exists( T& var_to_set, const std::string& attribute, T lower_bound ) { if (has_attribute(attribute)) { var_to_set = get_attribute(attribute, lower_bound); } } template void YAMLConfigHandle::set_if_attribute_exists( int& var_to_set, const std::string& attribute, T lower_bound ) { if (has_attribute(attribute)) { var_to_set = get_attribute(attribute, lower_bound); } } } //end namespace io } //end namespace specmicp #endif // SPECMICP_UTILS_SAFECONFIG_HPP diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index d687cfa..538d2e3 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,463 +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.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, true); + 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 swapping ({ - {"Si(OH)4", "SiO(OH)3[-]"} + {"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(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.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-9); + 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, true); + 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-8)); + 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 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)); } }