diff --git a/CMakeLists.txt b/CMakeLists.txt index 69b6804..a7a4440 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,173 +1,176 @@ project(specmicp) cmake_minimum_required(VERSION 2.8) ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package find_package(CxxTest) # Necessary if you want the tests ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### set(PROJECT_TEST_DIR ${PROJECT_SOURCE_DIR}/tests) set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${EIGEN3_INCLUDE_DIR}) include_directories(${PROJECT_SOURCE_DIR}) ################ 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}/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) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/thermodata.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # -------- Database ---------------- set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${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 ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/data_container.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.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 ----------- # 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) # --------- Test ---------------------- # 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) 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) diff --git a/tests/specmicp/alkaliactivated.cpp b/tests/specmicp/alkaliactivated.cpp new file mode 100644 index 0000000..32d0782 --- /dev/null +++ b/tests/specmicp/alkaliactivated.cpp @@ -0,0 +1,97 @@ + +#include +#include "specmicp/reaction_path.hpp" + +#include "utils/log.hpp" + +void solve_lot_reaction_path() +{ + specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; + + specmicp::database::Database database("data/cemdata_specmicp.js"); + std::shared_ptr data = database.get_database(); + std::map swapping ({ + {"H[+]","HO[-]"}, + {"Al[3+]","Al(OH)4[-]"}, + {"Si(OH)4", "SiO(OH)3[-]"} + }); + database.swap_components(swapping); + + std::shared_ptr model = std::make_shared(); + double mult = 2; + double m_naoh = 0.1; + double m_sio2 = mult*0.5; + double m_al2o3 = mult*0.5; + double m_cao = mult*1.0; + double m_gypsum = 0.2; + //double wc = 0.5; + double m_water = 1.0; + //double delta_h2co3 = 0.1; + double delta_cao = 0.1; + model->amount_aqueous = { + {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water -2*m_al2o3-m_cao-m_sio2, +delta_cao)}, + //{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)}, + {"HO[-]", specmicp::reaction_amount_t(m_naoh - 2*m_al2o3 +2*m_cao-m_sio2, -2*delta_cao)}, + {"Na[+]", specmicp::reaction_amount_t(m_naoh, 0)}, + {"Al(OH)4[-]", specmicp::reaction_amount_t(2*m_al2o3, 0)}, + {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_sio2, 0)}, + {"Ca[2+]", specmicp::reaction_amount_t(m_cao, -delta_cao)} + }; + model->amount_minerals = { + {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)} + }; + model->database_path = "data/cemdata_specmicp.js"; + + + double totiter = 0; + double totfact = 0; + specmicp::ReactionPathDriver driver(model, data); + driver.dissolve_to_components(); + + Eigen::VectorXd totaq(data->nb_component); + Eigen::VectorXd x(data->nb_component+data->nb_mineral); + x(0) = 1.0; + x.block(1, 0, data->nb_component-1, 1).setConstant(-2); + x.block(data->nb_component, 0, data->nb_mineral-1, 1).setZero(); + + //driver.get_options().solver_options.penalization_factor =1.0; + driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; + //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; + driver.get_options().solver_options.use_scaling = false; + driver.get_options().solver_options.max_factorization_step = 2; + driver.get_options().solver_options.factor_gradient_search_direction = 50; + driver.get_options().solver_options.maxstep = 50; + driver.get_options().solver_options.maxiter_maxstep = 50; + driver.get_options().solver_options.max_iter = 100; + driver.get_options().allow_restart = true; + std::cout << "Reaction_path return_code nb_iter nb_fact pH"; + for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it) + { + std::cout << " " << *it; + } + for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it) + { + std::cout << " " << *it; + } + std::cout << std::endl; + int max_step=11; + for (int i=0; inb_component-1,1).transpose() << " " + << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << std::endl; + totiter += perf.nb_iterations; + totfact += perf.nb_factorization; + } + std::cout << "Average iterations : " << totiter/max_step << std::endl; + std::cout << "Average factorization : " << totfact/max_step << std::endl; +} + +int main() +{ + solve_lot_reaction_path(); +} diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp index b1e811f..1fd9696 100644 --- a/tests/specmicp/thermocarbo.cpp +++ b/tests/specmicp/thermocarbo.cpp @@ -1,162 +1,158 @@ #include #include "specmicp/reduced_system.hpp" #include "micpsolver/micpsolver.hpp" #include "micpsolver/micpsolver_min.hpp" #include "specmicp/reduced_system_solver.hpp" #include "database/database.hpp" #include "specmicp/reaction_path.hpp" void solve_lot_reaction_path() { specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; Eigen::VectorXd x(10); x(0) = 1.0; x(1) = -2.0; x(2) = -2.0; x(3) = -2.0; x(4) = -2.0; x(5) = 0.0; x(6) = 0.0; x(7) = 0.0; x(8) = 0.0; x(9) = 0.0; specmicp::database::Database database("data/cemdata_specmicp.js"); std::shared_ptr data = database.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); database.swap_components(swapping); std::shared_ptr model = std::make_shared(); double m_c3s = 0.7; double m_c2s = 0.3; double wc = 1.0; double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3; double delta_h2co3 = 0.1; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)}, {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)}, {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)} }; model->database_path = "data/cemdata_specmicp.js"; Eigen::VectorXd totaq(5); double totiter = 0; double totfact = 0; specmicp::ReactionPathDriver driver(model, data); driver.dissolve_to_components(); //driver.get_options().solver_options.penalization_factor =1.0; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 1; driver.get_options().solver_options.factor_gradient_search_direction = 50; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = true; std::cout << "Reaction_path return_code nb_iter nb_fact pH H2O HO C Ca Si " << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl; int max_step=41; for (int i=0; i data = database.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); database.swap_components(swapping); std::shared_ptr model = std::make_shared(); double m_c3s = 0.7; double m_c2s = 0.3; double wc = 1.0; double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3; double delta_h2co3 = 0.1; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)}, {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)}, {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)} }; //x(0) = m_water; model->database_path = "data/cemdata_specmicp.js"; - Eigen::VectorXd totaq(5); specmicp::ReactionPathDriver driver(model, data); driver.dissolve_to_components(); + + Eigen::VectorXd totaq(data->nb_component); + Eigen::VectorXd x(data->nb_component+data->nb_mineral); + x(0) = 1.0; + x.block(1, 0, data->nb_component-1, 1).setConstant(-2); + x.block(data->nb_component, 0, data->nb_mineral-1, 1).setZero(); + driver.get_options().solver_options.penalization_factor = 0.8; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 1; driver.get_options().solver_options.factor_gradient_search_direction = 50; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = false; driver.get_options().solver_options.power_descent_condition = 2.0; std::cout << "Reaction_path return_code nb_iter nb_fact pH H2O HO C Ca Si " << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl; specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); driver.total_aqueous_concentration(x, totaq); - std::cout << 0*(0.1) << " " << (int) perf.return_code << " " + std::cout << 0.0 << " " << (int) perf.return_code << " " << perf.nb_iterations << " " << perf.nb_factorization << " " << 14+x(1) << " " - << x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(5, 0, 5, 1).transpose() << std::endl; - + << x(0) << " " << totaq.block(1,0,data->nb_component-1,1).transpose() << " " + << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << std::endl; std::cout << "Solution \n" << x << std::endl; } int main() { specmicp::logger::ErrFile::stream() = &std::cerr; //solve(); solve_lot_reaction_path(); //for (int i=0; i<5000; ++i) solve_lot_reaction_path(); }