diff --git a/CMakeLists.txt b/CMakeLists.txt index aa717b8..89d98d7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,281 +1,291 @@ project(specmicp) cmake_minimum_required(VERSION 2.8.8) # Options : # - SPECMICP_USE_OPENMP : bool : use openmp to parallelize code # # ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Boost REQUIRED) include_directories(${Boost_INCLUDE_DIR}) find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package include_directories(${EIGEN3_INCLUDE_DIR}) if(EXISTS "${EIGEN3_INCLUDE_DIR}/unsupported/Eigen/src/IterativeSolvers/GMRES.h") add_definitions(-DEIGEN_UNSUPPORTED_FOUND) INCLUDE_DIRECTORIES("${EIGEN3_INCLUDE_DIR}/unsupported/") endif() ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11 -fuse-ld=gold") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") if(SPECMICP_USE_OPENMP) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") add_definitions(-DSPECMICP_USE_OPENMP) endif() message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### # 3rd party set(JSONCPP_DIR ${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) add_custom_target(3rdparty_header SOURCES ${JSONCPP_DIR}/json/json.h ${JSONCPP_DIR}/json/json-forwards.h ) include_directories(${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${PROJECT_SOURCE_DIR}) add_custom_target(common SOURCES ${PROJECT_SOURCE_DIR}/common.hpp ${PROJECT_SOURCE_DIR}/database.hpp ) ################ Modules ########################## # -------- MiCPSolver ---------------- set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(specmisolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolverold.hpp ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # ------- SpecMiCP --------------- set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) set(KINETICS_DIR ${SPECMICP_DIR}/kinetics) set(ADIMENSIONAL_DIR ${SPECMICP_DIR}/adimensional) +set(PROBLEM_SOLVER_DIR ${SPECMICP_DIR}/problem_solver) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/simulation_box.hpp ${SPECMICP_DIR}/reduced_system_solver_structs.hpp + + ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp + ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp + + ${PROBLEM_SOLVER_DIR}/formulation.hpp + ${KINETICS_DIR}/kinetic_model.hpp ${KINETICS_DIR}/kinetic_variables.hpp ${KINETICS_DIR}/kinetic_system_solver_structs.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp #${SPECMICP_DIR}/speciation_solver.cpp #${SPECMICP_DIR}/simple_composition.cpp ${ADIMENSIONAL_DIR}/adimensional_system.cpp + ${ADIMENSIONAL_DIR}/adimensional_system_solver.cpp + + ${PROBLEM_SOLVER_DIR}/dissolver.cpp # kinetic ${KINETICS_DIR}/kinetic_system.cpp ${KINETICS_DIR}/kinetic_system_solver.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # -------- Database ---------------- set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/data_container.cpp ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ${JSONCPP_DIR}/jsoncpp.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # --------- ODE int ------------------ set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ${ODEINT_DIR}/runge_kutta_step_structs.hpp ) # --------- DFPMSolver ----------------- set(DFPMSOLVER_DIR ${PROJECT_SOURCE_DIR}/dfpmsolver) add_custom_target(dfpmsolver_incl SOURCES ${DFPMSOLVER_DIR}/driver.hpp ${DFPMSOLVER_DIR}/driver.inl ${DFPMSOLVER_DIR}/driver_structs.hpp ${DFPMSOLVER_DIR}/dfpm_program.hpp # # Parabolic solver ${DFPMSOLVER_DIR}/parabolic_driver.hpp ${DFPMSOLVER_DIR}/parabolic_driver.inl ${DFPMSOLVER_DIR}/parabolic_program.hpp ${DFPMSOLVER_DIR}/parabolic_structs.hpp ) # --------- ReactMiCP ----------------- set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_DIR}/common.hpp # ${REACTMICP_DIR}/micpsolvers/parabolicprogram.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver_structs.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.inl # ${REACTMICP_DIR}/mesh.hpp ${REACTMICP_DIR}/meshes/mesh1dfwd.hpp ${REACTMICP_DIR}/meshes/mesh1d.hpp ${REACTMICP_DIR}/meshes/uniform_mesh1d.hpp ${REACTMICP_DIR}/meshes/axisymmetric_uniform_mesh1d.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.inl ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver_structs.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters_base.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/boundary_conditions.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_DIR}/systems/diffusion/diffusion.cpp ${REACTMICP_DIR}/systems/diffusion/diffusion_numbering.cpp ${REACTMICP_DIR}/systems/diffusion/diffusion_secondary.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_variables.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/variables.cpp ) add_library(reactmicp ${REACTMICP_LIBFILES}) # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp ${UTILS_DIR}/sparse_solver.hpp ${UTILS_DIR}/sparse_solver_structs.hpp ${UTILS_DIR}/options_handler.hpp ${UTILS_DIR}/perfs_handler.hpp ) # ------- Physics ---------------- set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # ------- Data ------------------- set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # ------ Documentation ----------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) endif(DOXYGEN_FOUND) # --------- Python interface ---------- add_subdirectory( cython ) # --------- Test ---------------------- add_subdirectory( tests ) # --------- Examples ------------------ add_subdirectory( examples ) diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp index 3ee52ea..b890c45 100644 --- a/src/specmicp/adimensional/adimensional_system.cpp +++ b/src/specmicp/adimensional/adimensional_system.cpp @@ -1,599 +1,595 @@ /*------------------------------------------------------- - Module : specmicp - File : adim_system.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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 #include #include "adimensional_system.hpp" #include "utils/log.hpp" #include "../equilibrium_data.hpp" #include "physics/constants.hpp" #include "physics/laws.hpp" +#include "adimensional_system_solution.hpp" + namespace specmicp { constexpr scalar_t log10 = std::log(10.0); AdimensionalSystem::AdimensionalSystem(RawDatabasePtr ptrdata, const Vector &totconc, bool conservation_water, index_t charge_keeper): m_data(ptrdata), m_tot_conc(totconc), m_secondary_conc(ptrdata->nb_aqueous), m_loggamma(ptrdata->nb_component+ptrdata->nb_aqueous), m_charge_keeper(charge_keeper) { number_eq(conservation_water); m_secondary_conc.setZero(); m_loggamma.setZero(); } AdimensionalSystem::AdimensionalSystem( RawDatabasePtr ptrdata, const Vector& totconc, - const EquilibriumState& previous_solution, + const AdimensionalSystemSolution& previous_solution, bool conservation_water, index_t charge_keeper ): m_data(ptrdata), m_tot_conc(totconc), - m_secondary_conc(previous_solution.molality_secondary()), - m_loggamma(previous_solution.activity_coefficient()), + m_secondary_conc(previous_solution.secondary_molalities), + m_loggamma(previous_solution.log_gamma), + m_ionic_strength(previous_solution.ionic_strength), m_charge_keeper(charge_keeper) { number_eq(conservation_water); } void AdimensionalSystem::number_eq(bool water_equation) { index_t neq = 0; m_ideq.reserve(m_data->nb_component+m_data->nb_mineral); if (water_equation) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } for (index_t i: m_data->range_aqueous_component()) { // Remove components that does not exist if (m_tot_conc(i) == 0 and i!= 1) { INFO << "Component " << m_data->labels_basis[i] << "has total concentration equal to zero, we desactivate it" << std::endl; m_nonactive_component.push_back(i); m_ideq.push_back(no_equation); continue; } else { m_ideq.push_back(neq); ++neq; } } m_nb_free_vars = neq; for (index_t m: m_data->range_mineral()) { bool can_precipitate = true; // Remove minerals that cannot precipitate for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_data->nu_mineral(m, *it) != 0) { can_precipitate = false; break; // this is not a mineral that can precipitate } } if (can_precipitate) { m_ideq.push_back(neq); ++neq; } else { m_ideq.push_back(no_equation); } } m_nb_tot_vars = neq; m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars; // aqueous species m_active_aqueous.reserve(m_data->nb_aqueous); for (index_t j: m_data->range_aqueous()) { bool can_exist = true; for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it) { if (m_data->nu_aqueous(j,*it) != 0) { can_exist = false; } } m_active_aqueous.push_back(can_exist); } } scalar_t AdimensionalSystem::saturation_water(const Vector& x) const { if (ideq_w() != no_equation) return x(ideq_w()); else return 1.0 - sum_saturation_minerals(x); } scalar_t AdimensionalSystem::saturation_mineral(const Vector& x, index_t mineral) const { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral); if (ideq_min(mineral) == no_equation) return 0.0; else return x(ideq_min(mineral)); } scalar_t AdimensionalSystem::sum_saturation_minerals(const Vector& x) const { scalar_t sum_saturations = 0.0; for (index_t mineral: m_data->range_mineral()) { sum_saturations += saturation_mineral(x, mineral); } return sum_saturations; } scalar_t AdimensionalSystem::density_water() const { scalar_t scaling = 1.0; if (mass_unit() == units::MassUnit::gram) scaling *= 1e-3; if (length_unit() == units::LengthUnit::centimeter) scaling *= 1e-6; return scaling*1e3; } scalar_t AdimensionalSystem::molar_volume_mineral(index_t mineral) const { return m_data->molar_volume_mineral(mineral, length_unit()); } scalar_t AdimensionalSystem::residual_water(const Vector& x) const { const scalar_t conc_w = density_water()*saturation_water(x); scalar_t res = m_tot_conc(0) - conc_w/m_data->molar_mass_basis_si(0); for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) continue; res -= conc_w*m_data->nu_aqueous(j, 0)*m_secondary_conc(j); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation or m_data->nu_mineral(m, 0) == 0.0) continue; res -= m_data->nu_mineral(m, 0)*saturation_mineral(x, m)/molar_volume_mineral(m); } return res; } double AdimensionalSystem::residual_component(const Vector &x, index_t component) const { const scalar_t conc_w = density_water()*saturation_water(x); scalar_t res = m_tot_conc(component) - conc_w*component_molality(x, component); for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) continue; res -= conc_w*m_data->nu_aqueous(j, component)*secondary_molality(j); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; res -= m_data->nu_mineral(m, component)*saturation_mineral(x, m)/molar_volume_mineral(m); } return res; } double AdimensionalSystem::residual_mineral(const Vector& x, index_t m) const { scalar_t res = m_data->logk_mineral(m); for (index_t i: m_data->range_aqueous_component()) { if (m_data->nu_mineral(m, i) != 0) res -= m_data->nu_mineral(m, i)*(log_component_molality(x, i) + log_gamma_component(i)); } return res; } double AdimensionalSystem::residual_charge_conservation(const Vector& x) const { scalar_t res = 0.0; for (index_t i: m_data->range_aqueous_component()) { if (m_data->charge_component(i) != 0 and ideq_paq(i) != no_equation) res += m_data->charge_component(i)*component_molality(x, i); } for (index_t j: m_data->range_aqueous()) { if (m_data->charge_aqueous(j) == 0 and not is_aqueous_active(j)) continue; res += m_data->charge_aqueous(j)*secondary_molality(j); } return res; } void AdimensionalSystem::get_residuals(const Vector& x, Vector& residual) { set_secondary_concentration(x); if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x); for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) continue; if (is_charge_keeper(i)) residual(ideq_paq(i)) = residual_charge_conservation(x); else residual(ideq_paq(i)) = residual_component(x, i); } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m); } } // //non-optimized Finite difference, for test only //void AdimensionalSystem::get_jacobian(Eigen::VectorXd& x, // Eigen::MatrixXd& jacobian) //{ // const int neq = total_variables(); // Eigen::VectorXd res(total_variables()); // Eigen::VectorXd perturbed_res(total_variables()); // get_residuals(x, res); // for (int j=0; jmolar_mass_basis_si(0); for (index_t j: m_data->range_aqueous()) { if ( m_data->nu_aqueous(j, 0) == 0 and not is_aqueous_active(j)) continue; tmp -= m_data->nu_aqueous(j, 0)*secondary_molality(j); } jacobian(0, 0) = tmp; const scalar_t conc_w = density_water()*saturation_water(x); for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp = 0; for (index_t j: m_data->range_aqueous()) { tmp -= m_data->nu_aqueous(j,0)*m_data->nu_aqueous(j,k)*secondary_molality(j); } jacobian(0, ideq_paq(k)) = conc_w*log10 * tmp; } for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(0, ideq_min(m)) = -m_data->nu_mineral(m, 0)/molar_volume_mineral(m); } } } void AdimensionalSystem::jacobian_aqueous_components(Vector& x, Matrix& jacobian) { for (index_t i: m_data->range_aqueous_component()) { const index_t idp = ideq_paq(i); if (idp == no_equation) continue; // Charge balance equation // ---------------------- if (is_charge_keeper(i)) { // Aqueous components for (index_t k: m_data->range_aqueous_component()) { const index_t idc = ideq_paq(k); if (idc == no_equation) continue; scalar_t tmp_drdb = 0.0; if (m_data->charge_component(k) != 0.0) tmp_drdb = m_data->charge_component(k); // Secondary species for (index_t j: m_data->range_aqueous()) { if ( not is_aqueous_active(j) or m_data->nu_aqueous(j, k) == 0.0 or m_data->charge_aqueous(j) == 0.0 ) continue; scalar_t tmp_value = m_data->nu_aqueous(j, k)*m_data->charge_aqueous(j); tmp_value *= secondary_molality(j)/component_molality(x, k); tmp_drdb += tmp_value; } jacobian(idp, idc) += component_molality(x,k)*log10*tmp_drdb; } } // Mass balance equation // --------------------- else { const scalar_t conc_w = density_water()*saturation_water(x); // Aqueous components for (index_t k: m_data->range_aqueous_component()) { if (ideq_paq(k) == no_equation) continue; scalar_t tmp_iip = 0; if (k == i) tmp_iip -= component_molality(x, i)*log10; for (index_t j: m_data->range_aqueous()) { tmp_iip -= log10*m_data->nu_aqueous(j, i)*m_data->nu_aqueous(j, k)*secondary_molality(j); } jacobian(idp, ideq_paq(k)) = conc_w*tmp_iip; } // Minerals for (index_t m: m_data->range_mineral()) { if (ideq_min(m) == no_equation) continue; jacobian(idp, ideq_min(m)) = - m_data->nu_mineral(m, i)/molar_volume_mineral(m); } // Water if (ideq_w() != no_equation) { scalar_t tmp_iw = -component_molality(x, i); for (index_t j: m_data->range_aqueous()) { if ( m_data->nu_aqueous(j, i) == 0 ) continue; tmp_iw -= m_data->nu_aqueous(j, i)*secondary_molality(j); } jacobian(idp, ideq_w()) = density_water()*tmp_iw; } } } } void AdimensionalSystem::jacobian_minerals(Vector& x, Matrix& jacobian) { for (index_t m: m_data->range_mineral()) { const index_t idm = ideq_min(m); if (idm == no_equation) continue; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) continue; jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i); } } } void AdimensionalSystem::set_secondary_concentration(const Vector& x) { for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) { m_secondary_conc(j) = 0; continue; } scalar_t logconc = -m_data->logk_aqueous(j) - m_loggamma(j+m_data->nb_component); for (index_t k: m_data->range_aqueous_component()) { if (m_data->nu_aqueous(j, k) == 0) continue; logconc += m_data->nu_aqueous(j, k)*(log_component_molality(x, k) + m_loggamma(k)); } m_secondary_conc(j) = pow10(logconc); } } void AdimensionalSystem::set_ionic_strength(const Vector& x) { scalar_t ionic = 0; for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation or m_data->charge_component(i) == 0) continue; ionic += component_molality(x, i)*std::pow(m_data->charge_component(i),2); } for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j) or m_data->charge_aqueous(j) == 0) continue; ionic += secondary_molality(j)*std::pow(m_data->charge_aqueous(j),2); } ionic_strength() = ionic/2; } void AdimensionalSystem::compute_log_gamma(const Vector& x) { set_ionic_strength(x); const scalar_t sqrti = std::sqrt(ionic_strength()); for (index_t i: m_data->range_aqueous_component()) { if (ideq_paq(i) == no_equation) { log_gamma_component(i) = 0.0; continue; } log_gamma_component(i) = laws::extended_debye_huckel( ionic_strength(), sqrti, m_data->charge_component(i), m_data->a_debye_component(i), m_data->b_debye_component(i) ); } for (index_t j: m_data->range_aqueous()) { if (not is_aqueous_active(j)) { log_gamma_secondary(j) = 0.0; continue; } log_gamma_secondary(j) = laws::extended_debye_huckel( ionic_strength(), sqrti, m_data->charge_aqueous(j), m_data->a_debye_aqueous(j), m_data->b_debye_aqueous(j) ); } } bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual) { not_in_linesearch = true; scalar_t previous_norm = m_loggamma.norm(); bool may_have_converged = false; if (norm_residual < nb_free_variables()) { // Use fixed point iterations for non-ideality for (int i=0; i<10; ++i) { set_secondary_concentration(x); compute_log_gamma(x); if (std::abs(previous_norm - m_loggamma.norm())/previous_norm < 1e-10) {may_have_converged = true; break;} previous_norm = m_loggamma.norm(); } } return may_have_converged; } double AdimensionalSystem::max_lambda(const Vector& x, const Vector& update) { if (ideq_w() != no_equation) { return 1.0/std::max(1.0, -update(0)/(0.9*x(0))); } else { return 1.0; } } -void AdimensionalSystem::reasonable_starting_guess(Vector &x, bool for_min) +void AdimensionalSystem::reasonable_starting_guess(Vector &x) { x.resize(m_data->nb_component+ m_data->nb_mineral); x(0) = 1.0; for (index_t i: m_data->range_aqueous_component()) { x(i) = -4.0; } - if (for_min) - x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero(); - else - x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero(); + x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero(); m_loggamma.setZero(); m_secondary_conc.setZero(); } void AdimensionalSystem::reasonable_restarting_guess(Vector& x) { x(0) = 1.0; for (index_t i=m_data->nb_component; inb_component+m_data->nb_mineral; ++i) { x(i) = 0.0; } m_loggamma.setZero(); m_secondary_conc.setZero(); } -EquilibriumState AdimensionalSystem::get_solution(const Vector& xtot, const Vector& x) +AdimensionalSystemSolution AdimensionalSystem::get_solution(const Vector& xtot, const Vector& x) { double previous_norm = m_loggamma.norm(); set_secondary_concentration(x); compute_log_gamma(x); if (std::abs(previous_norm - m_loggamma.norm()) > 1e-6) { WARNING << "Activity coefficient have not converged !" << std::endl << "output can not be trusted\n Difference : " +std::to_string(std::abs(previous_norm - m_loggamma.norm())); } - - return EquilibriumState(xtot, m_secondary_conc, - m_loggamma, m_ionic_strength, - m_data); - + return AdimensionalSystemSolution(xtot, m_secondary_conc, m_loggamma, m_ionic_strength); } void AdimensionalSystem::reduce_mineral_problem(Vector& true_var) { for (index_t mineral: m_data->range_mineral()) { if (ideq_min(mineral) == no_equation) continue; if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve) { if (true_var(ideq_min(mineral)) == 0.0) continue; for (index_t component: m_data->range_component()) { if (m_data->nu_mineral(mineral, component) == 0.0) continue; m_tot_conc(component) -= m_data->nu_mineral(mineral, component)*true_var(ideq_min(mineral)); } true_var(ideq_min(mineral)) = 0.0; } } } void AdimensionalSystem::reset_mineral_system(Vector& true_var, Vector& output_var) { for (index_t mineral: m_data->range_mineral()) { if (ideq_min(mineral) == no_equation) continue; if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve) { output_var(mineral) += true_var(ideq_min(mineral)); } } } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp index 8c16133..9c20bc3 100644 --- a/src/specmicp/adimensional/adimensional_system.hpp +++ b/src/specmicp/adimensional/adimensional_system.hpp @@ -1,318 +1,319 @@ /*------------------------------------------------------- - Module : specmicp - - File : adim_system.hpp + - File : adimensional_system.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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_ADIMSYSTEM_HPP -#define SPECMICP_ADIMSYSTEM_HPP +#ifndef SPECMICP_ADIMENSIONALSYSTEM_HPP +#define SPECMICP_ADIMENSIONALSYSTEM_HPP //! \file reduced_system.hpp The MiCP program to solve speciation #include "common.hpp" #include "database.hpp" #include "micpsolver/micpprog.hpp" #include "../simulation_box.hpp" #include "physics/units.hpp" + //! \namespace specmicp main namespace for the SpecMiCP solver namespace specmicp { -class EquilibriumState; +class AdimensionalSystemSolution; //! \brief The equilibrium speciation problem //! //! Represent the equilibrium speciation problem as a Mixed Complementarity program //! It is called reduced becaused it does not solve explicitely for the secondary species. //! class AdimensionalSystem : public micpsolver::MiCPProg, public units::UnitBaseClass { public: //! \brief Create a reduced system //! //! \param ptrthermo Shared pointer to the thermodynamic data //! \param totconc Vector containing the total concentrations (in mol/volume) for each components AdimensionalSystem( RawDatabasePtr ptrdata, const Vector& total_concentration, bool conservation_water=true, index_t charge_keeper=no_species ); AdimensionalSystem( RawDatabasePtr ptrdata, const Vector& totconc, - const EquilibriumState& previous_solution, + const AdimensionalSystemSolution& previous_solution, bool conservation_water=true, index_t charge_keeper=no_species ); // Variables // ========== //! \brief Return the total number of variables index_t total_variables() const {return m_nb_tot_vars;} //! \brief Return the number of 'free' variables (not subject to complementarity conditions) index_t nb_free_variables() const {return m_nb_free_vars;} //! \brief Return the number of variables subjected to the complementarity conditions index_t nb_complementarity_variables() const {return m_nb_compl_vars;} // The linear system // ================== //! \brief Return the residual for the water conservation equation scalar_t residual_water(const Vector& x) const; //! \brief Return the residual for the conservation equation of component (i) //! //! For (i==0), water use the method : 'residual_water' double residual_component(const Vector& x, index_t i) const; //! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$ double residual_mineral(const Vector& x, index_t m) const; //! \brief Charge conservation double residual_charge_conservation(const Vector& x) const; //! \brief Return the residuals void get_residuals(const Vector& x, Vector& residual); //! \brief Return the jacobian void get_jacobian(Vector& x, Matrix& jacobian); // Getters // ####### // Equation id access // ================== //! \brief Return the equation id index_t ideq(index_t i) const {return m_ideq[i];} //! \brief Return the equation number for conservation of water index_t ideq_w() const {return m_ideq[0];} //! \brief Return the equation number of component 'i' index_t ideq_paq(index_t i) const { specmicp_assert(i < m_data->nb_component); return m_ideq[i]; } //! \brief Return the equation number of mineral 'm' index_t ideq_min(index_t m) const { specmicp_assert(m < m_data->nb_mineral); return m_ideq[m+m_data->nb_component]; } // Variables // ========= // Primary // ------- //! \brief Return the saturation of water scalar_t saturation_water(const Vector& x) const; //! \brief Return the density of water scalar_t density_water() const; //! \brief Return the saturation (mol/volume) of a mineral scalar_t saturation_mineral(const Vector& x, index_t mineral) const; //! \brief Return the volume of minerals scalar_t molar_volume_mineral(index_t mineral) const; //! \brief Return the sum of saturation of the minerals scalar_t sum_saturation_minerals(const Vector& x) const; //! \brief Return the log_10 of component concentration scalar_t log_component_molality(const Vector& x, index_t component) const { specmicp_assert(ideq_paq(component) != no_equation and component < m_data->nb_component); return x(ideq_paq(component)); } //! \brief Return the concentration of 'component' scalar_t component_molality(const Vector& x, index_t component) const { return pow10(log_component_molality(x, component)); } // Secondary // --------- //! \brief Return the concentration of secondary specis 'aqueous' scalar_t secondary_molality(index_t aqueous) const { if (not m_active_aqueous[aqueous]) return 0.0; return m_secondary_conc(aqueous); } //! \brief Return true if 'aqueous' is active //! //! i.e. Return true if all its component are present in the system bool is_aqueous_active(index_t aqueous) const { return m_active_aqueous[aqueous]; } //! \brief Return log_10(γ_i) where i is a component scalar_t log_gamma_component(index_t component) const { return m_loggamma(component); } //! \brief Return log_10(γ_j) where j is a secondary aqueous species scalar_t log_gamma_secondary(index_t aqueous) const { return m_loggamma(m_data->nb_component + aqueous); } //! \brief Return the ionic strength scalar_t ionic_strength() const { return m_ionic_strength; } //! \brief Return the total concentration of 'component' scalar_t total_concentration(index_t component) const {return m_tot_conc(component);} // Equilibrium State // ================= //! \brief Return the equilibrium state of the system, the 'Solution' of the speciation problem - EquilibriumState get_solution(const Vector& xtot, const Vector& x); + AdimensionalSystemSolution get_solution(const Vector& xtot, const Vector& x); // Algorithm // ######### //! \brief Compute the ionic strength void set_ionic_strength(const Vector& x); //! \brief Compute the activity coefficients void compute_log_gamma(const Vector& x); // For the normal algorithm void set_secondary_concentration(const Vector& x); //! \brief Compute the activity model, called by the solver at the beginning of an iteration bool hook_start_iteration(const Vector &x, scalar_t norm_residual); //! \brief Return a scale factor to avoid negative mass during Newton's iteration double max_lambda(const Vector &x, const Vector &update); //! \brief Return the component that does not exist in the solution (inactive dof) const std::vector& get_non_active_component() {return m_nonactive_component;} //! \brief A reasonable (well... maybe...) starting guess - void reasonable_starting_guess(Vector& x, bool for_min=false); + void reasonable_starting_guess(Vector& x); //! \brief A reasonable (maybe...) restarting guess void reasonable_restarting_guess(Vector& x); //! \brief Reduce the mineral system if one cannot dissolve void reduce_mineral_problem(Vector& true_var); //! \brief Reset the mineral system to the true problem void reset_mineral_system(Vector& true_var, Vector& output_var); //! \brief set the charge keeper to 'component' //! //! The equation for the charge keeper is replaced by the charge balance //! //! Set component to 'no_species' to solve charge conservation implicitely. //! In this case, the total concentration must conserve the charge. void set_charge_keeper(index_t component) { specmicp_assert(component >= 1 and component < m_data->nb_component); m_charge_keeper = component; } //! \brief Return true if 'component' is thecharge keeper bool is_charge_keeper(index_t component) { return (component == m_charge_keeper); } // Private Members and attributes // ############################## private: // Jacobian // ######## void jacobian_water(Vector& x, Matrix& jacobian); void jacobian_aqueous_components(Vector& x, Matrix& jacobian); void jacobian_minerals(Vector& x, Matrix& jacobian); // Equation numbering // ################## void number_eq(bool water_equation); // Setter // ###### // main variables // -------------- // they are set by the solver // Secondary variables // ------------------- //! \brief Return a reference to log_10(γ_i) where i is a component scalar_t& log_gamma_component(index_t component) { return m_loggamma(component); } //! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species scalar_t& log_gamma_secondary(index_t aqueous) { return m_loggamma(m_data->nb_component + aqueous); } //! \brief Return a reference to the ionic strength scalar_t& ionic_strength() { return m_ionic_strength; } // Attributes // ########## RawDatabasePtr m_data; Vector m_tot_conc; index_t m_nb_tot_vars; index_t m_nb_free_vars; index_t m_nb_compl_vars; std::vector m_ideq; Vector m_secondary_conc; Vector m_gas_pressure; - scalar_t m_ionic_strength; Vector m_loggamma; + scalar_t m_ionic_strength; index_t m_charge_keeper; bool not_in_linesearch; std::vector m_nonactive_component; std::vector m_active_aqueous; std::vector m_active_gas; }; } // end namespace specmicp -#endif // SPECMICP_REDUCEDSYSTEM_HPP +#endif // SPECMICP_ADIMENSIONALSYSTEM_HPP diff --git a/src/specmicp/adimensional/adimensional_system_solution.hpp b/src/specmicp/adimensional/adimensional_system_solution.hpp new file mode 100644 index 0000000..d5fd664 --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solution.hpp @@ -0,0 +1,34 @@ +#ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP +#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP + +#include "common.hpp" + +namespace specmicp { + +//! \brief Solution of an adimensional system +struct AdimensionalSystemSolution +{ + Vector main_variables; //!< The main variables + Vector secondary_molalities; //!< Molalities of secondary aqueous + Vector log_gamma; //!< Log_10 of activities coefficients + scalar_t ionic_strength; //!< Ionic strength + + AdimensionalSystemSolution() {} + + AdimensionalSystemSolution( + const Vector& variables, + const Vector& aqueous_molalities, + const Vector& loggama, + scalar_t the_ionic_strength + ): + main_variables(variables), + secondary_molalities(aqueous_molalities), + log_gamma(loggama), + ionic_strength(the_ionic_strength) + {} + +}; + +} // end namespace specmicp + +#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp new file mode 100644 index 0000000..4aebad7 --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solver.cpp @@ -0,0 +1,215 @@ +/*------------------------------------------------------- + + - Module : specmicp + - File : reduced_system_solver + - Author : Fabien Georget + +Copyright (c) 2014, Fabien 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: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * 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. + * Neither the name of the Princeton University 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 OWNER 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_solver.hpp" +#include "micpsolver/micpsolver.hpp" + +#include "specmicp/adimensional/adimensional_system_solution.hpp" + +namespace specmicp { + + +AdimensionalSystemSolver::AdimensionalSystemSolver( + RawDatabasePtr data, + const Vector& totconc, + AdimensionalSystemSolverOptions options + ): + OptionsHandler(options), + m_data(data), + m_system(std::make_shared( + data, + totconc, + options.conservation_water, + options.charge_keeper + )), + m_var(Vector::Zero(data->nb_component+data->nb_mineral)) +{} + +AdimensionalSystemSolver::AdimensionalSystemSolver(RawDatabasePtr data, + const Vector& totconc, + const AdimensionalSystemSolution& previous_solution, + AdimensionalSystemSolverOptions options + ): + OptionsHandler(options), + m_data(data), + m_system(std::make_shared( + data, + totconc, + previous_solution, + options.conservation_water, + options.charge_keeper + )), + m_var(Vector::Zero(data->nb_component+data->nb_mineral)) +{} + +AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution(const Eigen::VectorXd& x) +{ + set_true_variable_vector(x); + return m_system->get_solution(x, m_var); +} + + +micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init) +{ + m_system->set_units(get_options().units_set); + if (init) + { + m_system->reasonable_starting_guess(x); + } + set_true_variable_vector(x); + micpsolver::MiCPPerformance perf = solve(); + if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized and get_options().allow_restart) + { + WARNING << "Failed to solve the system ! - We shake it up and start again"; + const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor; + const scalar_t save_factor_descent = get_options().solver_options.factor_descent_condition; + get_options().solver_options.factor_descent_condition *= 1e-2; + if (save_penalization_factor == 1) + get_options().solver_options.penalization_factor = 0.8; + set_return_vector(x); + m_system->reasonable_restarting_guess(x); + set_true_variable_vector(x); + micpsolver::MiCPPerformance perf2 = solve(); + get_options().solver_options.penalization_factor = save_penalization_factor; + get_options().solver_options.factor_descent_condition = save_factor_descent; + perf += perf2; + } + if (perf.return_code > micpsolver::MiCPSolverReturnCode::NotConvergedYet) set_return_vector(x); + return perf; +} + +void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x) +{ + const std::vector& non_active_component = m_system->get_non_active_component(); + if (non_active_component.size() == 0) + { + // we still copy the data, if we failed to solve the problem, we can restart + if (get_options().conservation_water) + m_var = x; // direct copy + else + m_var = x.block(1, 0, x.rows()-1, 1); + for (int i=0; irange_aqueous_component()) + { + auto it = std::find(non_active_component.begin(), non_active_component.end(),i); + if (it != non_active_component.end()) continue; + m_var(new_i) = x(i); + ++new_i; + } + for (index_t m: m_data->range_mineral()) + { + for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it) + { + if (m_data->nu_mineral(m, *it) != 0) continue; + m_var(new_i) = x(m_data->nb_component+m); + ++new_i; + } + } + m_var.conservativeResize(new_i); + } + m_system->reduce_mineral_problem(m_var); +} + +micpsolver::MiCPPerformance AdimensionalSystemSolver::solve() +{ + + micpsolver::MiCPSolver solver(m_system); + solver.set_options(get_options().solver_options); + solver.solve(m_var); + return solver.get_performance(); + +} + +void AdimensionalSystemSolver::set_return_vector(Vector& x) +{ + const std::vector& non_active_component = m_system->get_non_active_component(); + if (non_active_component.size() == 0) + { + if (get_options().conservation_water) + x = m_var; //direct copy + else + x.block(1, 0, x.rows()-1, 1) = m_var; + // at that point we should have the correct solution + } + else + { + uindex_t new_i = 0; + if (get_options().conservation_water) + { + x(0) = m_var(new_i); + ++new_i; + } + for (index_t i: m_data->range_aqueous_component()) + { + auto it = std::find(non_active_component.begin(), non_active_component.end(),i); + if (it != non_active_component.end()) + { + x(i) = -HUGE_VAL; + continue; + } + x(i) = m_var(new_i) ; + ++new_i; + } + for (index_t m: m_data->range_mineral()) + { + for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it) + { + if (m_data->nu_mineral(m, *it) != 0) + { + x(m_data->nb_component+m) = 0; + continue; + } + x(m_data->nb_component+m) =m_var(new_i); + ++new_i; + } + } + } + m_system->reset_mineral_system(m_var, x); +} + +} // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solver.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp new file mode 100644 index 0000000..9612aae --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solver.hpp @@ -0,0 +1,111 @@ +/*------------------------------------------------------- + + - Module : specmicp + - File : reduced_system_solver.hpp + - Author : Fabien Georget + +Copyright (c) 2014, Fabien 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: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * 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. + * Neither the name of the Princeton University 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 OWNER 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_ADIMENSIONALSYSTEMSOLVER_HPP +#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP + +#include "database.hpp" +#include "adimensional_system.hpp" +#include "adimensional_system_solver_structs.hpp" +#include "utils/options_handler.hpp" + +//! \file reduced_system_solver.hpp Solve a reduced system + +namespace specmicp { + +// forward declaration +class AdimensionalSystemSolution; + +//! \brief Solver a reduced system +//! +//! Take care of possibly non-existing component in the system +class AdimensionalSystemSolver: public OptionsHandler +{ +public: + using SystemPtr = std::shared_ptr; + + //! Default constructor - you need to use another method to initialize it correctly + AdimensionalSystemSolver() {} + + AdimensionalSystemSolver(RawDatabasePtr data, + const Vector& totconc, + AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions() + ); + + AdimensionalSystemSolver(RawDatabasePtr data, + const Vector& totconc, + const AdimensionalSystemSolution& previous_solution, + AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions() + ); + + //! \brief solve the problem using initial guess x + //! + //! \param[in,out] x in -> initial guess, out -> solution + //! \param init if true, the algorithm guess a starting point + micpsolver::MiCPPerformance solve(Vector& x, bool init=false); + + //! \brief Return the system used for the computation + SystemPtr get_system() {return m_system;} + + //! \brief Return the solution in a manageable form + AdimensionalSystemSolution get_raw_solution(const Eigen::VectorXd& x); + +private: + //! \brief set up the true variable vector + void set_true_variable_vector(const Vector& x); + //! \brief set up the true solution vector + //! + //! add zero components + void set_return_vector(Vector& x); + //! \brief solve the problem + micpsolver::MiCPPerformance solve(); + + RawDatabasePtr m_data; //! The raw database + SystemPtr m_system; //! The system to solve + Vector m_var; //! Copy of the solution vector (necessary in case of failing) +}; + +//! \brief Solve a reduced system, function provided for convenience +inline micpsolver::MiCPPerformance solve_reduced_system( + std::shared_ptr data, + Vector& totconc, + Vector& x + ) +{ + return AdimensionalSystemSolver(data, totconc).solve(x); +} + + +} // end namespace specmicp + +#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP diff --git a/src/specmicp/adimensional/adimensional_system_solver_structs.hpp b/src/specmicp/adimensional/adimensional_system_solver_structs.hpp new file mode 100644 index 0000000..47767f8 --- /dev/null +++ b/src/specmicp/adimensional/adimensional_system_solver_structs.hpp @@ -0,0 +1,33 @@ +#ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP +#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP + +#include "common.hpp" +#include "micpsolver/micpsolver_structs.hpp" +#include "physics/units.hpp" + +namespace specmicp { + +//! \brief Options for the Equilibrium solver +//! +//! Most of the options are contained in the MiCP solver options +struct AdimensionalSystemSolverOptions +{ + micpsolver::MiCPSolverOptions solver_options; //!< Options of the MiCP solver + bool conservation_water; //!< Solve the conservation of liquid water + index_t charge_keeper; //!< Species that explicitely conserve the charge + bool allow_restart; //!< Allow the restarting if the problem failed the first part + units::UnitsSet units_set; //!< Set of units + + AdimensionalSystemSolverOptions(): + solver_options(), + conservation_water(true), + charge_keeper(no_species), + allow_restart(true) + { + } +}; + + +} // end namespace specmicp + +#endif //SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP diff --git a/src/specmicp/problem_solver/dissolver.cpp b/src/specmicp/problem_solver/dissolver.cpp new file mode 100644 index 0000000..40b8668 --- /dev/null +++ b/src/specmicp/problem_solver/dissolver.cpp @@ -0,0 +1,130 @@ +#include "dissolver.hpp" + +#include "formulation.hpp" + +namespace specmicp { + +Vector Dissolver::dissolve(const Formulation& the_problem) +{ + for (index_t component: m_database->range_component()) + { + m_components_to_keep[component] = false; + } + m_total_concentration.setZero(); + + set_solvent(the_problem.mass_solution); + for (auto it=the_problem.concentration_aqueous.begin(); it!=the_problem.concentration_aqueous.end(); ++it) + { + dissolve_aqueous(it->first, it->second, the_problem.mass_solution); + } + for (auto it=the_problem.amount_minerals.begin(); it!=the_problem.amount_minerals.end(); ++it) + { + dissolve_mineral(it->first, it->second); + } + keep_extra_components(the_problem.extra_components_to_keep); + reduce_problem(); + return m_total_concentration; +} + + +void Dissolver::set_solvent(scalar_t amount) +{ + index_t idw = m_database_hl.component_label_to_id("H2O"); + specmicp_assert(idw != no_species and idw == 0); + m_total_concentration(0) += amount/m_database->molar_mass_basis(idw, mass_unit()); + m_components_to_keep[0] = true; +} + +void Dissolver::dissolve_aqueous(std::string label, scalar_t concentration, scalar_t mass_solution) +{ + index_t idaq = m_database_hl.component_label_to_id(label); + if (idaq != no_species) + { + m_total_concentration(idaq) += mass_solution*concentration; + m_components_to_keep[idaq] = true; + } + else + { + index_t idsaq = m_database_hl.aqueous_label_to_id(label); + if (idsaq == no_species) + { + throw std::invalid_argument("Unknown species : " + label + "."); + } + for (index_t component: m_database->range_component()) + { + if (m_database->nu_aqueous(idsaq, component) == 0.0) continue; + m_total_concentration(component) += + m_database->nu_aqueous(idsaq, component)*mass_solution*concentration; + m_components_to_keep[component] = true; + } + } +} + +void Dissolver::dissolve_mineral(std::string label, scalar_t amount) +{ + index_t id_mineral = m_database_hl.mineral_label_to_id(label); + if (id_mineral != no_species) + { + for (index_t component: m_database->range_component()) + { + if (m_database->nu_mineral(id_mineral, component) == 0.0) continue; + m_total_concentration(component) += + m_database->nu_mineral(id_mineral, component)*amount; + m_components_to_keep[component] = true; + } + } + else + { + index_t id_minkin = m_database_hl.mineral_kinetic_label_to_id(label); + if (id_minkin == no_species) + { + throw std::invalid_argument("Unknown mineral : " + label + "."); + } + for (index_t component: m_database->range_component()) + { + if (m_database->nu_mineral_kinetic(id_minkin, component) == 0.0) continue; + m_total_concentration(component) += + m_database->nu_mineral_kinetic(id_minkin, component)*amount; + m_components_to_keep[component] = true; + } + } +} + +void Dissolver::keep_extra_components(const std::vector& list_component_to_keep) +{ + for (auto it=list_component_to_keep.begin(); it != list_component_to_keep.end(); ++it) + { + index_t idc = m_database_hl.component_label_to_id(*it); + if (idc == no_species) + { + throw std::invalid_argument("Species '" + *it + "' is not a component"); + } + m_components_to_keep[idc] = true; + } +} + +void Dissolver::reduce_problem() +{ + std::vector to_remove; + to_remove.reserve(m_database->nb_component); + index_t new_id = 0; + Vector reduced_tot_conc(m_database->nb_component); + for (index_t component: m_database->range_component()) + { + if (m_components_to_keep[component]) + { + reduced_tot_conc(new_id) = m_total_concentration(component); + ++new_id; + } + else + { + to_remove.push_back(component); + } + } + m_database_hl.remove_components(to_remove); + reduced_tot_conc.conservativeResize(m_database->nb_component); + m_total_concentration = reduced_tot_conc; + +} + +} // end namespace specmicp diff --git a/src/specmicp/problem_solver/dissolver.hpp b/src/specmicp/problem_solver/dissolver.hpp new file mode 100644 index 0000000..4ea990f --- /dev/null +++ b/src/specmicp/problem_solver/dissolver.hpp @@ -0,0 +1,44 @@ +#ifndef SPECMICP_SPECMICP_PROBLEMSOLVER_DISSOLVER_HPP +#define SPECMICP_SPECMICP_PROBLEMSOLVER_DISSOLVER_HPP + +#include "common.hpp" +#include "database.hpp" +#include "physics/units.hpp" + +namespace specmicp { + +class Formulation; + +class Dissolver: public units::UnitBaseClass +{ +public: + + Dissolver(RawDatabasePtr the_database): + m_database(the_database), + m_database_hl(m_database), + m_total_concentration(the_database->nb_component), + m_components_to_keep(the_database->nb_component) + { + m_total_concentration.setZero(); + } + + Vector dissolve(const Formulation& the_problem); + +private: + void set_solvent(scalar_t amount); + void dissolve_aqueous(std::string label, scalar_t concentration, scalar_t mass_solution); + void dissolve_mineral(std::string label, scalar_t amount); + void keep_extra_components(const std::vector& list_component_to_keep); + void reduce_problem(); + + + RawDatabasePtr m_database; + Database m_database_hl; + Vector m_total_concentration; + std::vector m_components_to_keep; + +}; + +} // end namespace specmicp + +#endif //SPECMICP_SPECMICP_PROBLEMSOLVER_DISSOLVER_HPP diff --git a/src/specmicp/problem_solver/formulation.hpp b/src/specmicp/problem_solver/formulation.hpp new file mode 100644 index 0000000..b79bb47 --- /dev/null +++ b/src/specmicp/problem_solver/formulation.hpp @@ -0,0 +1,22 @@ +#ifndef SPECMICP_SPECMICP_PROBLEMSOLVER_FORMULATION_HPP +#define SPECMICP_SPECMICP_PROBLEMSOLVER_FORMULATION_HPP + +#include "common.hpp" +#include +#include +#include + +namespace specmicp { + +struct Formulation +{ + scalar_t mass_solution; + std::map concentration_aqueous; + std::map amount_minerals; + + std::vector extra_components_to_keep; +}; + +} // end namespace specmicp + +#endif // SPECMICP_SPECMICP_PROBLEMSOLVER_FORMULATION_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e1b8391..2ce628b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,164 +1,165 @@ ##################### External packages ######################### find_package(CxxTest) # Necessary if you want the tests # catch for unit testing include(ExternalProject) find_package(Git REQUIRED) ExternalProject_Add( catch PREFIX ${CMAKE_BINARY_DIR}/catch GIT_REPOSITORY https://github.com/philsquared/Catch.git TIMEOUT 10 UPDATE_COMMAND ${GIT_EXECUTABLE} pull CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" LOG_DOWNLOAD ON ) ################## Configuration ################################ set(PROJECT_TEST_DIR ${CMAKE_CURRENT_SOURCE_DIR}) ExternalProject_Get_Property(catch source_dir) set(CATCH_INCLUDE_DIR ${source_dir}/include CACHE INTERNAL "Path to include folder for Catch") include_directories(${CATCH_INCLUDE_DIR}) include_directories(${PROJECT_TEST_DIR}) ################## Tests ######################################## enable_testing(true) # CXX Test # ======== if(CXXTEST_FOUND) set(CXXTEST_USE_PYTHON TRUE) include_directories(${CXXTEST_INCLUDE_DIR}) enable_testing() # MiCP Solver # ------------ CXXTEST_ADD_TEST(test_cond_number test_cond_number.cpp ${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp) CXXTEST_ADD_TEST(test_ncp_funtion test_ncp_function.cpp ${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp) CXXTEST_ADD_TEST(test_micpsolver test_micpsolver.cpp ${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp) # Database # -------- CXXTEST_ADD_TEST(test_data_reader test_data_reader.cpp ${PROJECT_TEST_DIR}/database/test_data_reader.hpp) target_link_libraries(test_data_reader specmicp_database jsoncpp) CXXTEST_ADD_TEST(test_data_selector test_data_selector.cpp ${PROJECT_TEST_DIR}/database/test_data_selector.hpp) target_link_libraries(test_data_selector specmicp_database jsoncpp) # Odeint # ------ CXXTEST_ADD_TEST(test_butchertableau test_butchertableau.cpp ${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp) CXXTEST_ADD_TEST(test_embeddedrungekutta test_embeddedrungekutta.cpp ${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp) endif() # Test 'Perso' # ============ add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp) target_link_libraries(thermocarbo specmicp specmicp_database jsoncpp) add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp) target_link_libraries(carboalu specmicp specmicp_database jsoncpp) add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp) target_link_libraries(alkaliactivated specmicp specmicp_database jsoncpp) add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp) target_link_libraries(rate_nicoleau specmicp specmicp_database jsoncpp) add_executable(kinetics ${PROJECT_TEST_DIR}/specmicp/kinetics.cpp) target_link_libraries(kinetics specmicp specmicp_database jsoncpp) set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv ) add_custom_target(data_nicoleau SOURCES ${DATA_TEST}) file(INSTALL ${DATA_TEST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/) # Catch test #=========== # Global step # ----------- add_custom_target(catch_header SOURCES test_database.hpp) add_executable(test_reactmicp_diffusion ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/test_reactmicp_diffusion.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/diffusion.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/numbering.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/secondary.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/solving.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/utils.cpp ) target_link_libraries(test_reactmicp_diffusion reactmicp specmicp specmicp_database) # Mesh # ---- add_executable(test_meshes ${PROJECT_TEST_DIR}/reactmicp/meshes/test_meshes.cpp ${PROJECT_TEST_DIR}/reactmicp/meshes/test_uniform_mesh_1d.cpp ${PROJECT_TEST_DIR}/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp ) # Saturated Diffusion # ------------------- add_executable(test_reactmicp_saturated_diffusion ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/test_reactmicp_saturated_diffusion.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/variables.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/solver.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/transport_program.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/reactive_transport.cpp ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/utils.cpp ) target_link_libraries(test_reactmicp_saturated_diffusion reactmicp specmicp specmicp_database) # Specmicp : Adim system # ---------------------- add_executable(test_specmicp_adim ${PROJECT_TEST_DIR}/specmicp/adim/test_specmicp_adim.cpp ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system.cpp ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_solver.cpp + ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_problem_solver.cpp ) target_link_libraries(test_specmicp_adim specmicp specmicp_database) diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp new file mode 100644 index 0000000..3bd6802 --- /dev/null +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -0,0 +1,98 @@ +#include "catch.hpp" + +#include +#include "utils/log.hpp" + +#include "specmicp/adimensional/adimensional_system_solver.hpp" +#include "specmicp/adimensional/adimensional_system_solution.hpp" +#include "specmicp/problem_solver/formulation.hpp" +#include "specmicp/problem_solver/dissolver.hpp" + + +TEST_CASE("Solving adimensinal system with problem setting", "[specmicp, MiCP, program, adimensional, solver]") { + + specmicp::logger::ErrFile::stream() = &std::cerr; + specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; + + + SECTION("Solving problem with dissolver - simpler problem") { + + specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); + std::map swapping ({ + {"H[+]","HO[-]"}, + {"Si(OH)4", "SiO(OH)3[-]"} + }); + thedatabase.swap_components(swapping); + specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); + + specmicp::Formulation formulation; + formulation.mass_solution = 0.8; + formulation.amount_minerals = {{"Portlandite", 10.0}}; + + specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); + total_concentrations *= 1e3; + + specmicp::AdimensionalSystemSolver solver(raw_data, total_concentrations); + + solver.get_options().charge_keeper = 1; + solver.get_options().solver_options.maxstep = 50.0; + solver.get_options().solver_options.maxiter_maxstep = 100; + solver.get_options().solver_options.use_crashing = false; + solver.get_options().solver_options.use_scaling = false; + solver.get_options().solver_options.factor_descent_condition = 1e-10; + solver.get_options().solver_options.factor_gradient_search_direction = 200; + + specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral); + x.setConstant(0.0); + x(0) = 0.8; + x.segment(1, raw_data->nb_component-1).setConstant(-2); + + specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true); + + REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); + specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + + REQUIRE(x(0) == Approx(0.8).epsilon(1e-4)); + REQUIRE(x(3) == Approx(0.32947).epsilon(1e-2)); + + } + + SECTION("Solving problem with dissolver") { + + specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); + std::map swapping ({ + {"H[+]","HO[-]"}, + {"Si(OH)4", "SiO(OH)3[-]"} + }); + thedatabase.swap_components(swapping); + specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); + + specmicp::Formulation formulation; + formulation.mass_solution = 0.8; + formulation.amount_minerals = {{"C3S", 0.7}, {"C2S", 0.3}}; + formulation.extra_components_to_keep = {"HCO3[-]", }; + + specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); + total_concentrations *= 1e3; + total_concentrations(2) = 500; + + specmicp::AdimensionalSystemSolver solver(raw_data, total_concentrations); + + solver.get_options().charge_keeper = 1; + solver.get_options().solver_options.maxstep = 50.0; + solver.get_options().solver_options.maxiter_maxstep = 100; + solver.get_options().solver_options.use_crashing = false; + solver.get_options().solver_options.use_scaling = false; + solver.get_options().solver_options.factor_descent_condition = 1e-10; + solver.get_options().solver_options.factor_gradient_search_direction = 200; + + specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral); + + specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true); + + REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); + + specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); + + } +} diff --git a/tests/specmicp/adim/adimensional_system_solver.cpp b/tests/specmicp/adim/adimensional_system_solver.cpp index cacf3ea..2de125e 100644 --- a/tests/specmicp/adim/adimensional_system_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_solver.cpp @@ -1,108 +1,144 @@ #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" + specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); return thedatabase.get_database(); } TEST_CASE("Solving adimensinal system", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; SECTION("Solving simple case") { specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); specmicp::Vector total_concentrations(thedatabase->nb_component); total_concentrations << 3.0e4, 2e4 , 1e4; specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 100; options.maxiter_maxstep = 100; options.use_crashing = false; options.use_scaling = true; options.penalization_factor = 0.8; options.max_factorization_step = 4; options.factor_descent_condition = 1e-10; std::shared_ptr system= std::make_shared(thedatabase, total_concentrations); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral); x << 0.8, -2.0, -2.3, 0.0; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); REQUIRE(x(0) == Approx(0.54045).epsilon(1e-4)); REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volum in cm^3") { specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); specmicp::Vector total_concentrations(thedatabase->nb_component); total_concentrations << 0.03, 0.02, 0.01; specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 100; options.maxiter_maxstep = 100; options.use_crashing = false; - options.use_scaling = true; + options.use_scaling = false; options.penalization_factor = 0.8; options.max_factorization_step = 4; options.factor_descent_condition = 1e-10; std::shared_ptr system= std::make_shared(thedatabase, total_concentrations); system->length_unit() = specmicp::units::LengthUnit::centimeter; specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral); x << 0.8, -2.0, -2.3, 0.0; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); REQUIRE(x(0) == Approx(0.54045).epsilon(1e-4)); REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4)); } + + SECTION("Automatic solver") { + + specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); + + specmicp::Vector total_concentrations(thedatabase->nb_component); + total_concentrations << 0.03, 0.02, 0.01; + specmicp::Vector x(thedatabase->nb_component+thedatabase->nb_mineral); + x << 0.8, -2.0, -2.3, 0.0; + + + specmicp::AdimensionalSystemSolver solver(thedatabase, total_concentrations); + + solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; + solver.get_options().solver_options.maxstep = 100.0; + solver.get_options().solver_options.maxiter_maxstep = 100; + solver.get_options().solver_options.use_crashing = false; + solver.get_options().solver_options.use_scaling = false; + solver.get_options().solver_options.factor_descent_condition = 1e-10; + solver.get_options().solver_options.factor_gradient_search_direction = 100; + + solver.solve(x); + + + REQUIRE(x(0) == Approx(0.54045).epsilon(1e-4)); + REQUIRE(x(1) == Approx(-1.48059).epsilon(1e-4)); + REQUIRE(x(2) == Approx(-1.8538).epsilon(1e-4)); + REQUIRE(x(3) == Approx(0.32966).epsilon(1e-4)); + + } + }