diff --git a/src/specmicp/CMakeLists.txt b/src/specmicp/CMakeLists.txt index e648f42..44d89ec 100644 --- a/src/specmicp/CMakeLists.txt +++ b/src/specmicp/CMakeLists.txt @@ -1,110 +1,113 @@ # The speciation solver # ===================== # Directories set(ADIMENSIONAL_DIR adimensional) set(ADIMKINETICS_DIR adimensional_kinetics) set(PROBLEM_SOLVER_DIR problem_solver) # Headers file without a source file add_custom_target(specmicp_incl SOURCES ${ADIMENSIONAL_DIR}/adimensional_system_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_numbering.hpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm_structs.hpp ${PROBLEM_SOLVER_DIR}/formulation.hpp ${ADIMKINETICS_DIR}/kinetic_model.hpp ${ADIMKINETICS_DIR}/kinetic_variables.hpp ${ADIMKINETICS_DIR}/kinetic_system_solver_structs.hpp ) # Source for the library set(SPECMICP_LIBFILE ${ADIMENSIONAL_DIR}/adimensional_system.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solver.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solution_extractor.cpp + ${ADIMENSIONAL_DIR}/adimensional_system_solution_saver.cpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm.cpp ${ADIMENSIONAL_DIR}/equilibrium_curve.cpp ${PROBLEM_SOLVER_DIR}/dissolver.cpp # kinetic ${ADIMKINETICS_DIR}/kinetic_system.cpp ${ADIMKINETICS_DIR}/kinetic_system_solver.cpp ${ADIMKINETICS_DIR}/kinetic_system_euler_solver.cpp + + ${JSONCPP_DIR}/jsoncpp.cpp ) # -fvisibility=hidden set_visibility_hidden_flag( ${SPECMICP_LIBFILE}) if(SPECMICP_DEBUG_EQUATION_FD_JACOBIAN) set_source_files_properties(${ADIMENSIONAL_DIR}/adimensional_system.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_DEBUG_EQUATION_FD_JACOBIAN" ) endif() add_library(specmicp SHARED ${SPECMICP_LIBFILE}) install(TARGETS specmicp LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # include files # -------------- set(SPECMICP_ADIMENSIONAL_INCLUDES ${ADIMENSIONAL_DIR}/adimensional_system_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_numbering.hpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution_extractor.hpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm.hpp ${ADIMENSIONAL_DIR}/equilibrium_curve.hpp ) set(SPECMICP_PROBLEM_SOLVER_INCLUDES ${PROBLEM_SOLVER_DIR}/dissolver.hpp ${PROBLEM_SOLVER_DIR}/formulation.hpp ) set(SPECMICP_ADIMENSIONAL_KINETICS_INCLUDES ${ADIMKINETICS_DIR}/kinetic_model.hpp ${ADIMKINETICS_DIR}/kinetic_variables.hpp ${ADIMKINETICS_DIR}/kinetic_system_solver_structs.hpp ${ADIMKINETICS_DIR}/kinetic_system.hpp ${ADIMKINETICS_DIR}/kinetic_system_solver.hpp ${ADIMKINETICS_DIR}/kinetic_system_euler_solver.hpp ) install(FILES ${SPECMICP_ADIMENSIONAL_INCLUDES} DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp/adimensional ) install(FILES ${SPECMICP_PROBLEM_SOLVER_INCLUDES} DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp/problem_solver ) install(FILES ${SPECMICP_ADIMENSIONAL_KINETICS_INCLUDES} DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp/adimensional_kinetics ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(specmicp_static STATIC ${SPECMICP_LIBFILE}) install(TARGETS specmicp_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(specmicp_static EXCLUDE_FROM_ALL STATIC ${SPECMICP_LIBFILE}) endif() set_target_properties(specmicp_static PROPERTIES OUTPUT_NAME specmicp) diff --git a/src/specmicp/adimensional/adimensional_system_solution_saver.cpp b/src/specmicp/adimensional/adimensional_system_solution_saver.cpp new file mode 100644 index 0000000..3630383 --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solution_saver.cpp @@ -0,0 +1,355 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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 "adimensional_system_solution_saver.hpp" + +#include "adimensional_system_solution.hpp" +#include "../../../3rdparty/jsoncpp/json/json.h" + +#include + +#define SECTION_MAIN "main_variables" +#define SECTION_LOGGAMMA "log_gamma" +#define SECTION_GAS "gas_fugacities" +#define SECTION_SORBED "sorbed_molalities" + +#define VALUE_FREE_SURFACE "free_surface" +#define VALUE_IONIC_STRENGTH "ionic_strength" +#define VALUE_INERT "inert_vol_fraction" +#define VALUE_COMPONENT "components" +#define VALUE_AQUEOUS "aqueous" +#define VALUE_MINERAL "minerals" + +#include + +namespace specmicp { + +void check_for_mandatory_value(const Json::Value& root, + const std::string& attribute) +{ + if (not root.isMember(attribute)) + { + throw std::invalid_argument( + "Missing required attribute in solution : '" + attribute + "'."); + } +} + + +//! \brief Read a solution from a JSON file +AdimensionalSystemSolution parse_solution_json( + const std::string& filename, + const RawDatabasePtr& the_database + ) +{ + return AdimensionalSystemSolutionSaver(the_database).parse_solution(filename); +} + + +//! \brief Read a solution from a JSON file +AdimensionalSystemSolution parse_solution_json( + std::istream& input, + const RawDatabasePtr& the_database + ) +{ + return AdimensionalSystemSolutionSaver(the_database).parse_solution(input); +} + + +void AdimensionalSystemSolutionSaver::save_solution( + const std::string& filename, + const AdimensionalSystemSolution& solution + ) +{ + std::ofstream ofile(filename); + ofile << format_solution(solution); + ofile.close(); +} + +std::string AdimensionalSystemSolutionSaver::format_solution( + const AdimensionalSystemSolution& solution + ) +{ + Json::Value root; + root[SECTION_MAIN] = format_main_variables(solution); + root[SECTION_LOGGAMMA] = format_loggamma(solution); + root[VALUE_IONIC_STRENGTH] = solution.ionic_strength; + root[SECTION_GAS] = format_gas_fugacities(solution); + root[SECTION_SORBED] = format_sorbed_molalities(solution); + root[VALUE_INERT] = solution.inert_volume_fraction; + Json::StyledWriter writer; + return writer.write(root); + +} + + +AdimensionalSystemSolution +AdimensionalSystemSolutionSaver::parse_solution(std::istream& input) +{ + Json::Value root; + Json::Reader reader; + bool parsingSuccessful = reader.parse(input, root); + if ( !parsingSuccessful ) + { + // report to the user the failure and their locations in the document. + std::string message("Failed to parse configuration\n" + reader.getFormattedErrorMessages()); + throw std::runtime_error(message); + } + + + AdimensionalSystemSolution the_solution; + + read_main_variables(root, the_solution); + read_loggamma(root, the_solution); + read_gas_fugacities(root, the_solution); + read_sorbed_molalities(root, the_solution); + + check_for_mandatory_value(root, VALUE_INERT); + the_solution.inert_volume_fraction = root[VALUE_INERT].asDouble(); + + check_for_mandatory_value(root, VALUE_IONIC_STRENGTH); + the_solution.ionic_strength = root[VALUE_IONIC_STRENGTH].asDouble(); + + the_solution.is_valid = true; + return the_solution; +} + + +AdimensionalSystemSolution +AdimensionalSystemSolutionSaver::parse_solution(const std::string& filename) +{ + std::ifstream ofile(filename); + AdimensionalSystemSolution the_solution = parse_solution(ofile); + ofile.close(); + return the_solution; +} + + +Json::Value AdimensionalSystemSolutionSaver::format_main_variables( + const AdimensionalSystemSolution& solution) +{ + Json::Value main; + Json::Value components; + for (auto component: m_data->range_component()) + { + components[m_data->get_label_component(component)] = + solution.main_variables(dof_component(component)); + } + Json::Value minerals; + for (auto mineral: m_data->range_mineral()) + { + minerals[m_data->get_label_mineral(mineral)] = + solution.main_variables(dof_mineral(mineral)); + } + main[VALUE_COMPONENT] = std::move(components); + main[VALUE_FREE_SURFACE] = solution.main_variables(dof_surface()); + main[VALUE_MINERAL] = std::move(minerals); + return main; +} + +Json::Value AdimensionalSystemSolutionSaver::format_loggamma( + const AdimensionalSystemSolution& solution + ) +{ + Json::Value gamma; + Json::Value components; + for (auto component: m_data->range_aqueous_component()) + { + components[m_data->get_label_component(component)] = + solution.log_gamma(dof_component_gamma(component)); + } + Json::Value aqueous_list; + for (auto aqueous: m_data->range_aqueous()) + { + aqueous_list[m_data->get_label_aqueous(aqueous)] = + solution.log_gamma(dof_aqueous_gamma(aqueous)); + } + gamma[VALUE_COMPONENT] = std::move(components); + gamma[VALUE_AQUEOUS] = std::move(aqueous_list); + return gamma; +} + +Json::Value AdimensionalSystemSolutionSaver::format_gas_fugacities( + const AdimensionalSystemSolution& solution + ) +{ + Json::Value gas_fugacities; + for (auto gas: m_data->range_gas()) + { + gas_fugacities[m_data->get_label_gas(gas)] = + solution.gas_fugacities(gas); + } + return gas_fugacities; +} + + +Json::Value AdimensionalSystemSolutionSaver::format_sorbed_molalities( + const AdimensionalSystemSolution& solution + ) +{ + Json::Value sorbed_molalities; + for (auto sorbed: m_data->range_sorbed()) + { + sorbed_molalities[m_data->get_label_sorbed(sorbed)] = + solution.sorbed_molalities(sorbed); + } + return sorbed_molalities; +} + +void AdimensionalSystemSolutionSaver::read_main_variables( + const Json::Value& root, + AdimensionalSystemSolution& solution + ) +{ + // Initialization + // -------------- + solution.main_variables = Vector(total_dofs()); + solution.main_variables(dof_water()) = 0; + solution.main_variables(dof_electron()) = - INFINITY; + for (auto id: m_data->range_aqueous_component()) + { + solution.main_variables(dof_component(id)) = - INFINITY; + } + solution.main_variables(dof_surface()) = - INFINITY; + for (auto mineral: m_data->range_mineral()) + { + solution.main_variables(dof_mineral(mineral)) = 0; + } + // Parsing values + // -------------- + check_for_mandatory_value(root, SECTION_MAIN); + const Json::Value& main = root[SECTION_MAIN]; + // components + check_for_mandatory_value(main, VALUE_COMPONENT); + const Json::Value& components = main[VALUE_COMPONENT]; + for (auto& it: components.getMemberNames()) + { + index_t id = m_data->get_id_component(it); + if (id == no_species) + { + throw std::invalid_argument("Component '" + it + "' does not exist in the database !"); + } + solution.main_variables(dof_component(id)) = components[it].asDouble(); + } + // sorption + if (main.isMember(VALUE_FREE_SURFACE) and not main[VALUE_FREE_SURFACE].isNull()) + { + solution.main_variables(dof_surface()) = main[VALUE_FREE_SURFACE].asDouble(); + } + // minerals + check_for_mandatory_value(main, VALUE_MINERAL); + const Json::Value& minerals = main[VALUE_MINERAL]; + for (auto& it: minerals.getMemberNames()) + { + index_t id = m_data->get_id_mineral(it); + if (id == no_species) + { + throw std::invalid_argument("Mineral '" + it + "' does not exist in the database !"); + } + solution.main_variables(dof_mineral(id)) = minerals[it].asDouble(); + } +} + +void AdimensionalSystemSolutionSaver::read_loggamma( + const Json::Value& root, + AdimensionalSystemSolution& solution + ) +{ + solution.log_gamma = Vector::Zero(m_data->nb_component()+m_data->nb_aqueous()); + check_for_mandatory_value(root, SECTION_LOGGAMMA); + const Json::Value& loggamma = root[SECTION_LOGGAMMA]; + // components + check_for_mandatory_value(loggamma, VALUE_COMPONENT); + const Json::Value& components = loggamma[VALUE_COMPONENT]; + for (auto& it: components.getMemberNames()) + { + index_t id = m_data->get_id_component(it); + if (id == no_species) + { + throw std::invalid_argument("Component '" + it + "' does not exist in the database !"); + } + solution.log_gamma(dof_component_gamma(id)) = components[it].asDouble(); + } + // aqueous + check_for_mandatory_value(loggamma, VALUE_AQUEOUS); + const Json::Value& aqueous = loggamma[VALUE_AQUEOUS]; + for (auto& it: aqueous.getMemberNames()) + { + index_t id = m_data->get_id_aqueous(it); + if (id == no_species) + { + throw std::invalid_argument("Aqueous '" + it + "' does not exist in the database !"); + } + solution.log_gamma(dof_aqueous_gamma(id)) = aqueous[it].asDouble(); + } +} + +void AdimensionalSystemSolutionSaver::read_gas_fugacities( + const Json::Value& root, + AdimensionalSystemSolution& solution + ) +{ + solution.gas_fugacities = Vector::Zero(m_data->nb_gas()); + check_for_mandatory_value(root, SECTION_GAS); + if (root[SECTION_GAS].isNull()) return; + const Json::Value& gas_section = root[SECTION_GAS]; + for (auto& it: gas_section.getMemberNames()) + { + index_t id = m_data->get_id_gas(it); + if (id == no_species) + { + throw std::invalid_argument("Gas '" + it + "' does not exist in the database !"); + } + solution.gas_fugacities(id) = gas_section[it].asDouble(); + } +} + +void AdimensionalSystemSolutionSaver::read_sorbed_molalities( + const Json::Value& root, + AdimensionalSystemSolution& solution + ) +{ + solution.sorbed_molalities = Vector::Zero(m_data->nb_sorbed()); + check_for_mandatory_value(root, SECTION_SORBED); + if (root[SECTION_SORBED].isNull()) return; + const Json::Value& sorbed_section = root[SECTION_SORBED]; + for (auto& it: sorbed_section.getMemberNames()) + { + index_t id = m_data->get_id_sorbed(it); + if (id == no_species) + { + throw std::invalid_argument("sorbed species '" + it + "' does not exist in the database !"); + } + solution.sorbed_molalities(id) = sorbed_section[it].asDouble(); + } +} + +} //end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solution_saver.hpp b/src/specmicp/adimensional/adimensional_system_solution_saver.hpp new file mode 100644 index 0000000..d979171 --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solution_saver.hpp @@ -0,0 +1,143 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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_ADIMENSIONALSYSTEMSOLUTIONSAVER_HPP +#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONSAVER_HPP + +#include "../../types.hpp" +#include "../../database.hpp" +#include "../../../3rdparty/jsoncpp/json/json-forwards.h" + +#include "adimensional_system_numbering.hpp" + +#include + +namespace specmicp { + +struct AdimensionalSystemSolution; + +//! \brief This is a serializer for the AdimensionalSystemSolution +//! +//! This class format the solution is a JSON-formatted string. +class SPECMICP_DLL_PUBLIC AdimensionalSystemSolutionSaver : AdimemsionalSystemNumbering +{ +public: + AdimensionalSystemSolutionSaver(const RawDatabasePtr& the_database): + AdimemsionalSystemNumbering(the_database) + {} + + //! \brief Format the solution as a JSON-formatted string + std::string format_solution(const AdimensionalSystemSolution& solution); + + //! \brief Save the solution as a JSON file + void save_solution( + const std::string& filename, + const AdimensionalSystemSolution& solution + ); + + //! \brief Create a solution from a JSON input + AdimensionalSystemSolution parse_solution(std::istream& input); + + //! \brief Create a solution from a JSON file + AdimensionalSystemSolution parse_solution(const std::string& filename); + +private: + //! \brief Format the main variables + Json::Value SPECMICP_DLL_LOCAL format_main_variables( + const AdimensionalSystemSolution& solution); + + //! \brief Format the log of activity coefficients + Json::Value SPECMICP_DLL_LOCAL format_loggamma( + const AdimensionalSystemSolution& solution + ); + //! \brief Format the gas fugacities + Json::Value SPECMICP_DLL_LOCAL format_gas_fugacities( + const AdimensionalSystemSolution& solution + ); + //! \brief Format the sorbed molalities + Json::Value SPECMICP_DLL_LOCAL format_sorbed_molalities( + const AdimensionalSystemSolution& solution + ); + + void SPECMICP_DLL_LOCAL read_main_variables( + const Json::Value& root, + AdimensionalSystemSolution& solution + ); + void SPECMICP_DLL_LOCAL read_loggamma( + const Json::Value& root, + AdimensionalSystemSolution& solution + ); + void SPECMICP_DLL_LOCAL read_gas_fugacities( + const Json::Value& root, + AdimensionalSystemSolution& solution + ); + + void SPECMICP_DLL_LOCAL read_sorbed_molalities( + const Json::Value& root, + AdimensionalSystemSolution& solution + ); +}; + +//! \brief Serialize a solution as a JSON-formatted string +inline std::string format_solution_json( + const RawDatabasePtr& the_database, + const AdimensionalSystemSolution& the_solution + ) +{ + return AdimensionalSystemSolutionSaver(the_database).format_solution(the_solution); +} + +//! \brief Serialize a solution and save it in a file +inline void save_solution_json( + const std::string& filename, + const RawDatabasePtr& the_database, + const AdimensionalSystemSolution& the_solution + ) +{ + AdimensionalSystemSolutionSaver(the_database).save_solution(filename, the_solution); +} + +//! \brief Read a solution from a JSON file +AdimensionalSystemSolution parse_solution_json( + const std::string& filename, + const RawDatabasePtr& the_database + ); + +//! \brief Read a solution from a JSON file +AdimensionalSystemSolution parse_solution_json( + std::istream& input, + const RawDatabasePtr& the_database + ); + +} //end namespace specmicp + +#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONSAVER_HPP diff --git a/tests/specmicp/adim/adimensional_system_solver.cpp b/tests/specmicp/adim/adimensional_system_solver.cpp index 11508c4..0880b8f 100644 --- a/tests/specmicp/adim/adimensional_system_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_solver.cpp @@ -1,162 +1,176 @@ #include "catch.hpp" #include #include "utils/log.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "micpsolver/micpsolver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" +#include "specmicp/adimensional/adimensional_system_solution_saver.hpp" + #include "database/database.hpp" specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("Solving adimensional system", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Solving simple case") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 3.0e4; total_concentration(id_oh) = 2.0e4; total_concentration(id_ca) = 1.0e4; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 100; options.maxiter_maxstep = 100; options.disable_crashing(); options.disable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.48059).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.8538).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volum in cm^3") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 10; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); units::UnitsSet unit_set; unit_set.length = units::LengthUnit::centimeter; system->set_units(unit_set); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.48059).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.8538).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Automatic solver") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::Vector x; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints); solver.initialise_variables(x, 0.8, -2.0); //x(solver.dof_surface()) = -HUGE_VAL; solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; solver.get_options().solver_options.maxstep = 10.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.disable_descent_direction(); solver.get_options().solver_options.factor_gradient_search_direction = 100; solver.solve(x); AdimensionalSystemSolution solution = solver.get_raw_solution(x); AdimensionalSystemSolutionExtractor extr(solution, thedatabase, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4)); CHECK(extr.log_molality_component(id_oh) == Approx(-1.48059).epsilon(1e-4)); CHECK(extr.log_molality_component(id_ca) == Approx(-1.8538).epsilon(1e-4)); //CHECK(extr.free_surface_concentration() == -HUGE_VAL); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4)); + + AdimensionalSystemSolutionSaver saver(thedatabase); + saver.save_solution("test_solution.js", solution); + + AdimensionalSystemSolution solution2 = saver.parse_solution("test_solution.js"); + REQUIRE(solution.main_variables(0) == Approx(solution2.main_variables(0))); + REQUIRE(solution.main_variables(id_oh) == Approx(solution2.main_variables(id_oh))); + REQUIRE(solution.main_variables(id_ca) == Approx(solution2.main_variables(id_ca))); + REQUIRE(extr.volume_fraction_mineral(id_ch) == Approx(solution2.main_variables(extr.dof_mineral(id_ch)))); + REQUIRE(solution.log_gamma.norm() == Approx(solution2.log_gamma.norm())); + + } }