diff --git a/CMakeLists.txt b/CMakeLists.txt index ec612b5..1d015d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,246 +1,271 @@ #################################################### # # # SpecMiCP : CMakeLists.txt # # # #################################################### project(specmicp) cmake_minimum_required(VERSION 3.2) # CMake Options # ============= # For an explanation of the options see the INSTALL file option( SPECMICP_USE_OPENMP "Use OpenMP for parallelisation" ON ) option( SPECMICP_NO_DEBUG "Disable SpecMiCP assert" OFF ) option( SPECMICP_BUILD_STATIC "Build static libraries" OFF ) option( SPECMICP_BUILD_EXAMPLE "Build the examples" ON ) option( SPECMICP_BENCHMARK "Build benchmark" OFF ) option( SPECMICP_TEST "Enable testing" ON ) option( SPECMICP_BINARIES_USE_STATIC "Executables use static libraries" OFF ) # the following is only a debug options for developpers option( SPECMICP_DEBUG_EQUATION_FD_JACOBIAN "Use a finite difference jacobian" OFF ) # PGO sequence option( SPECMICP_PROFILE_GENERATE "Generate profile for PGO optimization" OFF ) option( SPECMICP_PROFILE_USE "Use profile for PGO optimization" OFF ) # LTO optimization option( SPECMICP_LTO "Use link time optimization" OFF ) option( SPECMICP_LD_GOLD "Use GNU gold linker" ON ) # global variables # ================ set(SPECMICP_VERSION 0.0.4) # External Package # ================ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake) # OpenMP #------- # not required but recommended if(SPECMICP_USE_OPENMP) find_package(OpenMP) + if (OPENMP_FOUND) + set(SPECMICP_HAVE_OPENMP ON) + endif(OPENMP_FOUND) endif() # Eigen # ----- find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package include_directories(${EIGEN3_INCLUDE_DIR}) # Eigen unsuported # GMRES.h is really the file we are using/looking for # If it doesn't exist then the solver will not be included in the list of the parse solvers set( EIGEN_GMRES_PATH "${EIGEN3_INCLUDE_DIR}/unsupported/Eigen/src/IterativeSolvers/GMRES.h") if( EXISTS ${EIGEN_GMRES_PATH} ) add_definitions(-DEIGEN_UNSUPPORTED_FOUND) include_directories("${EIGEN3_INCLUDE_DIR}/unsupported/") endif() # Yaml-cpp # -------- include(FindPkgConfig) pkg_check_modules(YAML REQUIRED yaml-cpp>=0.5) include_directories(${YAML_INCLUDE_DIRS}) link_directories(${YAML_LIBRARY_DIRS}) # HDF5 # ---- find_package(HDF5 REQUIRED COMPONENTS C CXX) # HDF5 is optional, must be checked by files that used it + +# glibc functions +# --------------- +include(CheckIncludeFile) +include(CheckFunctionExists) + +macro(check_required_include name var) + CHECK_INCLUDE_FILE( ${name} ${var} ) + if (NOT ${var}) + message(SEND_ERROR "Missing required include ${name}") + endif() +endmacro(check_required_include) + +check_required_include( "string.h" HAVE_STRING_H ) +check_required_include( "dirent.h" HAVE_DIRENT_H ) +check_required_include( "unistd.h" HAVE_UNISTD_H ) +check_required_include( "sys/stat.h" HAVE_SYS_STAT_H ) +check_required_include( "limits.h" HAVE_LIMITS_H ) +check_required_include( "stdlib.h" HAVE_STDLIB_H ) + +CHECK_FUNCTION_EXISTS( "secure_getenv" SPECMICP_HAVE_SECURE_GETENV ) + # sanitizer # --------- # to check memory, undefined behavior and all... include(SanitizerBuild) # compilation flags # ============================================================================ # check the availability of : # require C++11 standard by default set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) # -fuse-ld=gold include(gold_linker) # pgo optimization include(pg_optimization) # -fvisibility-hidden include(visibility_flag) # just a friendly warning if(NOT UNIX) message(WARNING "Not tested on non linux platform ! Probably won't work") endif() include(lto) # set the flags # ------------- if (OPENMP_FOUND) SET( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") endif() # let's be pedantic, it's always fun SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") # always best to check what we do... message(STATUS "c++ flags : ${CMAKE_CXX_FLAGS}") # Directories ######################################################################### # use gnu coding standards include(GNUInstallDirs) # the libraries install dir set( LIBRARY_INSTALL_DIR ${CMAKE_INSTALL_FULL_LIBDIR} CACHE PATH "Installation directory for libraries" ) # the static libraries install dir set( STATIC_LIBRARY_INSTALL_DIR ${CMAKE_INSTALL_FULL_LIBDIR} CACHE PATH "Installation directory for static libraries" ) # Binaries # -------- set( BIN_INSTALL_DIR ${CMAKE_INSTALL_FULL_BINDIR} CACHE PATH "Installation directory for the programs" ) # include #-------- set( INCLUDE_INSTALL_DIR ${CMAKE_INSTALL_FULL_INCLUDEDIR} CACHE PATH "Installation directory for the headers" ) # share #------ set( SHARE_INSTALL_DIR "${CMAKE_INSTALL_FULL_DATADIR}/specmicp/" CACHE PATH "Installation directory for the miscalleneous files..." ) mark_as_advanced( LIBRARY_INSTALL_DIR STATIC_LIBRARY_INSTALL_DIR BIN_INSTALL_DIR INCLUDE_INSTALL_DIR SHARE_INSTALL_DIR ) # CPP API : SpecMiCP / ReactMiCP ######################################################################### # the main source directory - the c++ api set(SPECMICP_CPP_API ${CMAKE_CURRENT_SOURCE_DIR}/src) include_directories( ${SPECMICP_CPP_API} ) add_subdirectory( ${SPECMICP_CPP_API} ) # the necessary libraries to link to set( SPECMICP_LIBS # just speciation solver specmicp specmicp_database specmicp_common ${YAML_LIBRARIES} ) set( REACTMICP_LIBS # reactive transport solver reactmicp dfpm ${SPECMICP_LIBS} ) # static versions set( SPECMICP_STATIC_LIBS specmicp_static specmicp_database_static specmicp_common_static ${YAML_LIBRARIES} ) set( REACTMICP_STATIC_LIBS reactmicp_static dfpm_static ${SPECMICP_STATIC_LIBS} ) # Databases ######################################################################### add_subdirectory( data ) # Documentation ######################################################################### # "common" documentation # ----------------------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) install(FILES ${DOCS_LIST} DESTINATION ${SHARE_INSTALL_DIR} ) # scripts # -------- add_subdirectory( scripts ) # Doxygen documentation # --------------------- add_subdirectory( doc ) # Tests ######################################################################### if( SPECMICP_TEST ) enable_testing(true) endif() add_subdirectory( tests ) # Examples ######################################################################### add_subdirectory( examples ) # Benchmark ######################################################################## if (SPECMICP_BENCHMARK) add_subdirectory( benchmark ) endif() # Docker ######################################################################## file(COPY docker/Dockerfile docker/run_tests.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/docker ) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b4290fb..148bb2f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,31 +1,22 @@ ###################### ## SPECMICP CPP API ## ###################### -# extra cmake macros +# configuration file +include_directories(${PROJECT_BINARY_DIR}/src/) -macro(src_use_openmp afile) - if(OPENMP_FOUND) - if(SPECMICP_USE_OPENMP) - set_property( SOURCE afile - APPEND PROPERTY - COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" - ) - endif(SPECMICP_USE_OPENMP) - endif(OPENMP_FOUND) -endmacro(src_use_openmp) #=============# # Libraries # #=============# add_subdirectory( specmicp_common ) add_subdirectory( specmicp_database ) add_subdirectory( specmicp ) add_subdirectory( dfpm ) add_subdirectory( reactmicp ) # Binaries ######################################################################### add_subdirectory( bin ) diff --git a/src/bin/reactmicp/CMakeLists.txt b/src/bin/reactmicp/CMakeLists.txt index 6b84fe1..c9d9422 100644 --- a/src/bin/reactmicp/CMakeLists.txt +++ b/src/bin/reactmicp/CMakeLists.txt @@ -1,32 +1,30 @@ # The ReactMiCP binaries set( UNSATURATED_BIN unsaturated.cpp ) -src_use_openmp(unsaturated.cpp) - add_executable(unsaturated_bin ${UNSATURATED_BIN}) set_target_properties(unsaturated_bin PROPERTIES OUTPUT_NAME reactmicp_unsaturated) if (UNSATURATED_BINARIES_USE_STATIC) target_link_libraries(unsaturated_bin reactmicp_static specmicp_static specmicp_database_static specmicp_common_static yaml-cpp ) else() target_link_libraries(unsaturated_bin reactmicp specmicp specmicp_database specmicp_common yaml-cpp ) endif() install(TARGETS unsaturated_bin DESTINATION ${BIN_INSTALL_DIR}) diff --git a/src/bin/reactmicp/unsaturated.cpp b/src/bin/reactmicp/unsaturated.cpp index 7bb31b0..0537590 100644 --- a/src/bin/reactmicp/unsaturated.cpp +++ b/src/bin/reactmicp/unsaturated.cpp @@ -1,702 +1,704 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ //! \file bin/reactmicp/unsaturated.cpp //! \brief Unsaturated #include "specmicp_common/types.hpp" #include "specmicp_common/cli/parser.hpp" #include "specmicp_common/plugins/plugin_manager.hpp" #include "specmicp_common/io/yaml.hpp" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/filesystem.hpp" #include "specmicp_common/string_algorithms.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp_database/io/configuration.hpp" #include "dfpm/io/configuration.hpp" #include "specmicp_common/io/configuration.hpp" #include "reactmicp/reactmicp_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "dfpm/solver/parabolic_structs.hpp" #include "reactmicp/io/hdf5_unsaturated.hpp" +#include "specmicp_common/config.h" + #include #include #include #include #include "specmicp_common/io/config_yaml_sections.h" #define PLUGIN_MODULE_UPSCALING "upscaling_factory" #define PLUGIN_MODULE_VARIABLES "variables_factory" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; const char* welcome = R"(=============================== = Reactmicp - unsaturated = =============================== (C) copyright 2014-2016 F.Georget, Princeton University version: not good enough yet )"; int main (int argc, char* argv[]); plugins::PluginManager& initialize_plugin_manager( io::YAMLConfigHandle& config ); SimulationData initialize_simul_data( io::YAMLConfigHandle& safe_config, const std::vector& db_dirs ); std::shared_ptr initialize_variables( std::shared_ptr raw_data, std::shared_ptr the_mesh, const units::UnitsSet& the_units, io::YAMLConfigHandle&& configuration ); struct Staggers { std::shared_ptr upscaling; std::shared_ptr chemistry; std::shared_ptr transport; }; Staggers initialize_staggers( SimulationData& simul_data, io::YAMLConfigHandle&& configuration ); std::shared_ptr initialize_upscaling_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::shared_ptr initialize_chemistry_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::shared_ptr initialize_transport_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::string initialize_output( io::YAMLConfigHandle&& conf ); int main(int argc, char *argv[]) { std::cout << welcome << std::endl; init_logger(&std::cout, logger::Warning); // CLI options cli::CommandLineParser parser; parser.add_option('i', "input", cli::ValueType::string, "input_file (yaml format)"); parser.add_option('d', "database_dir", std::string(""), "directory where database are stored"); parser.set_help_message("ReactMiCP Unsaturated System"); const auto ret_parse = parser.parse(argc, argv); if (ret_parse > 0) { return EXIT_SUCCESS; } std::string input_filepath = parser.get_option("input"); std::string database_dir = parser.get_option("database_dir"); std::vector db_dirs; const auto env = utils::get_env("SPECMICP_DATABASE_PATH"); if (not env.empty()) { db_dirs = utils::split(env, ';'); } if (not database_dir.empty()) {db_dirs.push_back(database_dir);} // Configuration File // ------------------ auto config = io::YAMLConfigFile::make(input_filepath); initialize_plugin_manager(*(config.get())); // Eigen initialization // -------------------- -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP Eigen::initParallel(); #endif // Simulation Data // --------------- SimulationData simul_data = initialize_simul_data(*config.get(), db_dirs); // Staggers // ------- Staggers staggers = initialize_staggers( simul_data, config->get_section(SPC_CF_S_STAGGERS) ); // Check the variables // ------------------- auto res = VariablesInterface(simul_data.variables).check_variables(); if (res != VariablesValidity::good) { const auto msg = "Variables were not initialized correctly." "Please see log for errors."; ERROR << msg; throw std::invalid_argument(msg); } // Solver // ------ solver::ReactiveTransportSolver react_solver( staggers.transport, staggers.chemistry, staggers.upscaling); io::configure_reactmicp_options( react_solver.get_options(), config->get_section(SPC_CF_S_REACTMICP) ); // Initialize output std::unique_ptr saver; solver::SimulationInformation simul_info = io::configure_simulation_information( config->get_section(SPC_CF_S_SIMULINFO)); solver::ReactiveTransportRunner runner(react_solver, 1.0, 10.0, // dummy value, set later simul_info); io::configure_reactmicp_timestepper( runner.get_timestepper_options(), config->get_section(SPC_CF_S_TIMESTEPPER) ); // if (config->has_section(SPC_CF_S_REACTOUTPUT)) { auto filepath = initialize_output( config->get_section(SPC_CF_S_REACTOUTPUT)); saver = make_unique( filepath, simul_data.variables, simul_data.units); auto* saver_ptr = saver.get(); auto out_pol = [saver_ptr]( scalar_t timestep, solver::VariablesBasePtr _) { saver_ptr->save_timestep(timestep); }; saver->save_timestep(0.0); runner.set_output_policy(out_pol); } // run info auto run_section = config->get_section(SPC_CF_S_RUN); auto run_until = run_section.get_required_attribute( SPC_CF_S_RUN_A_RUNUTNIL); // we destroy config to get information config.reset(nullptr); runner.run_until(run_until, simul_data.variables); return EXIT_SUCCESS; } // ==================== // // Simulation data // // ==================== // SimulationData initialize_simul_data( io::YAMLConfigHandle& safe_config, const std::vector& db_dirs ) { // Package the information SimulationData simul_data; // The units // --------- if (safe_config.has_section(SPC_CF_S_UNITS)) { simul_data.units = io::configure_units( safe_config.get_section(SPC_CF_S_UNITS)); } // The database // ------------ simul_data.raw_data = io::configure_database( safe_config.get_section(SPC_CF_S_DATABASE), db_dirs ); if (not simul_data.raw_data->is_valid()) { CRITICAL << "Invalid database from the configuration."; throw std::runtime_error("Database parsed from the configuration has" " been detected to be invalid. Abort."); } // The mesh // -------- simul_data.mesh1d = io::configure_mesh(safe_config.get_section(SPC_CF_S_MESH)); // The variables // ------------- simul_data.variables = initialize_variables( simul_data.raw_data, simul_data.mesh1d, simul_data.units, safe_config.get_section(SPC_CF_S_REACTMICPVARIABLES) ); if (not simul_data.raw_data->is_valid()) { CRITICAL << "Invalid database after variable initialization."; throw std::runtime_error("Database has been detected to be invalid" " after the variables initialization. Abort."); } // The database is supposed set after this point ! simul_data.raw_data->freeze_db(); // The boundary conditions // ----------------------- simul_data.bcs = io::configure_unsaturated_boundary_conditions( simul_data.mesh1d->nb_nodes(), simul_data.raw_data.get(), safe_config.get_section(SPC_CF_S_BOUNDARYCONDITIONS) ); // Simulation data ready return simul_data; } std::shared_ptr initialize_variables( std::shared_ptr raw_data, std::shared_ptr the_mesh, const units::UnitsSet& the_units, io::YAMLConfigHandle&& configuration ) { auto type = configuration.get_optional_attribute( SPC_CF_S_REACTMICPVARIABLES_A_TYPE, SPC_CF_V_PLUGIN); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN; if (type != SPC_CF_V_PLUGIN) { configuration.report_error( io::YAMLConfigError::InvalidArgument, "Only accepted value for attribute '" SPC_CF_S_REACTMICPVARIABLES_A_TYPE "' is '" SPC_CF_V_PLUGIN "'." ); } auto plugin_name = configuration.get_required_attribute( SPC_CF_A_PLUGIN_NAME ); plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); auto extra_plugin_file = configuration.get_optional_attribute( SPC_CF_A_PLUGIN_FILE, ""); if (extra_plugin_file != "") { auto retcode = plugin_manager.load_plugin(extra_plugin_file); if (not retcode) { CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl; std::runtime_error("Fail to load plugin file : " + extra_plugin_file + "."); } } auto initializer = plugin_manager.get_object( PLUGIN_MODULE_VARIABLES, plugin_name); if (initializer == nullptr) { CRITICAL << "No variable initializer found"; throw std::runtime_error("No variable initializer found"); } return initializer->initialize_variables( raw_data, the_mesh, the_units, configuration); } // ================ // // Plugin Manager // // ================ // plugins::PluginManager& initialize_plugin_manager( io::YAMLConfigHandle& config ) { plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); plugin_manager.register_module(PLUGIN_MODULE_VARIABLES, make_unique()); plugin_manager.register_module(PLUGIN_MODULE_UPSCALING, make_unique()); if (config.has_section(SPC_CF_S_PLUGINS)) { io::configure_plugin_manager(config.get_section(SPC_CF_S_PLUGINS)); } return plugin_manager; } // ================ // // Staggers // // ================ // Staggers initialize_staggers( SimulationData& simul_data, io::YAMLConfigHandle&& conf_staggers ) { Staggers staggers; auto* var_ptr = simul_data.variables.get(); staggers.upscaling = initialize_upscaling_stagger(simul_data, conf_staggers); staggers.upscaling->initialize(var_ptr); staggers.chemistry = initialize_chemistry_stagger(simul_data, conf_staggers); staggers.chemistry->initialize(var_ptr); staggers.transport = initialize_transport_stagger(simul_data, conf_staggers); staggers.transport->initialize(var_ptr); // check the database if (not simul_data.raw_data->is_valid()) { CRITICAL << "Change in database during stagger initialization !"; throw std::logic_error("Change in the database has been detected" " during the staggers initialization. Abort."); } return staggers; } std::shared_ptr initialize_upscaling_stagger( SimulationData& simul_data, io::YAMLConfigHandle& config_staggers ) { auto conf = config_staggers.get_section(SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER); auto type = conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE, SPC_CF_V_PLUGIN); // only plugin accepted for now if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN; if (type != SPC_CF_V_PLUGIN) { conf.report_error( io::YAMLConfigError::InvalidArgument, "Invalid argument for attribute '" SPC_CF_A_TYPE "' only value accepted is '" SPC_CF_V_PLUGIN "'."); } auto plugin_name = conf.get_required_attribute( SPC_CF_A_PLUGIN_NAME); plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); auto extra_plugin_file = conf.get_optional_attribute( SPC_CF_A_PLUGIN_FILE, ""); if (extra_plugin_file != "") { auto retcode = plugin_manager.load_plugin(extra_plugin_file); if (not retcode) { CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl; std::runtime_error("Fail to load plugin file : " + extra_plugin_file + "."); } } auto upscaling_factory = plugin_manager.get_object( PLUGIN_MODULE_UPSCALING, plugin_name); if (upscaling_factory == nullptr) { CRITICAL << "No upscaling stagger factory found"; throw std::runtime_error("No upscaling stagger factory found"); } return upscaling_factory->get_upscaling_stagger( simul_data, std::move(conf)); } std::shared_ptr initialize_chemistry_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ) { // no conf // ------- if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER)) { auto opts = EquilibriumOptionsVector::make( simul_data.mesh1d->nb_nodes()); return EquilibriumStagger::make( simul_data.variables, simul_data.bcs, opts ); } // conf is given // -------------- auto conf_chem = conf_staggers.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER); auto type = conf_chem.get_optional_attribute( SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE, SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM; if (type !=SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM) { CRITICAL << "Unknown chemistry stagger"; throw std::invalid_argument("Only equilibrium stagger supported."); } std::shared_ptr opts {nullptr}; if (conf_chem.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)) { opts = io::configure_unsaturated_equilibrium_options( simul_data.mesh1d->nb_nodes(), simul_data.units, conf_chem.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS) ); } else { opts = EquilibriumOptionsVector::make(simul_data.mesh1d->nb_nodes()); } return EquilibriumStagger::make( simul_data.variables, simul_data.bcs, opts ); } std::shared_ptr initialize_transport_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ) { const UnsaturatedVariables* const variables = simul_data.variables.get(); const database::DataContainer* const raw_data = simul_data.raw_data.get(); // No conf // ------- if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER)) { return UnsaturatedTransportStagger::make( simul_data.variables, simul_data.bcs, true, true); } // Conf is given // ------------- auto transport_conf = conf_staggers.get_section( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER); const bool merge_sat = transport_conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION, true); const bool merge_aq = transport_conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS, true); auto stagger = UnsaturatedTransportStagger::make( simul_data.variables, simul_data.bcs, merge_sat, merge_aq); // ! scaling => done in variables // default options // --------------- dfpmsolver::ParabolicDriverOptions default_opts; bool rewrite_default = false; if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS)) { io::configure_transport_options( default_opts, transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS) ); rewrite_default = true; } // Saturation option // ----------------- auto& sat_opts = *(stagger->get_saturation_options()); sat_opts = default_opts; if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS)) { io::configure_transport_options( sat_opts, transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS) ); } // Aqueous options // --------------- if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS)) { auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS); if (not conf.is_sequence()) { // config for every aqueous equations dfpmsolver::ParabolicDriverOptions opts; if (rewrite_default) opts = default_opts; io::configure_transport_options( opts, std::move(conf)); for (index_t aq: raw_data->range_aqueous_component()) { *(stagger->get_aqueous_options(aq)) = opts; } } else { // First set default is needed if (rewrite_default) { for (index_t aq: simul_data.raw_data->range_aqueous_component()) { *(stagger->get_aqueous_options(aq)) = default_opts; } } uindex_t size = conf.size(); for (uindex_t ind=0; ind( SPC_CF_A_COMPONENT); auto comp_id = raw_data->get_id_component(comp_label); io::configure_transport_options( *(stagger->get_aqueous_options(comp_id)), std::move(subconf) ); } } } // Gas options // ----------- if ( transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS) and not (merge_aq and merge_sat)) // skip conf if not needed { auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS); if (not conf.is_sequence()) { // config for every gas equations dfpmsolver::ParabolicDriverOptions opts; if (rewrite_default) opts = default_opts; io::configure_transport_options( opts, std::move(conf)); if ( not merge_sat and variables->component_has_gas(0)) { *(stagger->get_gas_options(0)) = opts; } if ( not merge_aq) { for (index_t aq: raw_data->range_aqueous_component()) { if (simul_data.variables->component_has_gas(aq) ) *(stagger->get_aqueous_options(aq)) = opts; } } } else { // First rewrite default if (rewrite_default) { // H2O if (not merge_sat and variables->component_has_gas(0)) { *(stagger->get_gas_options(0)) = default_opts; } // aqueous components if (not merge_aq) { for (index_t aq: raw_data->range_aqueous_component()) { if (simul_data.variables->component_has_gas(aq)) *(stagger->get_aqueous_options(aq)) = default_opts; } } } // Then go over all conf uindex_t size = conf.size(); for (uindex_t ind=0; ind( SPC_CF_A_COMPONENT); auto comp_id = raw_data->get_id_component(comp_label); // H2O if (comp_id == 0) { if (merge_sat) continue; io::configure_transport_options( *(stagger->get_gas_options(0)), std::move(subconf) ); } // aqueous component else { if (merge_aq) continue; io::configure_transport_options( *(stagger->get_gas_options(comp_id)), std::move(subconf) ); } } } } return stagger; } std::string initialize_output( io::YAMLConfigHandle&& conf ) { auto type = conf.get_optional_attribute( SPC_CF_S_REACTOUTPUT_A_TYPE, SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5; if (type != SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5) { conf.report_error( io::YAMLConfigError::InvalidArgument, "Only hdf5 is accepted for attribute " SPC_CF_S_REACTOUTPUT_A_TYPE " for now."); } auto filepath = conf.get_required_attribute( SPC_CF_S_REACTOUTPUT_A_FILEPATH); return filepath; } diff --git a/src/reactmicp/systems/saturated_react/CMakeLists.txt b/src/reactmicp/systems/saturated_react/CMakeLists.txt index d9031ef..3f7910a 100644 --- a/src/reactmicp/systems/saturated_react/CMakeLists.txt +++ b/src/reactmicp/systems/saturated_react/CMakeLists.txt @@ -1,33 +1,30 @@ set( reactmicp_saturated_srcs diffusive_upscaling_stagger.cpp equilibrium_stagger.cpp init_variables.cpp kinetic_stagger.cpp transport_program.cpp transport_stagger.cpp variables.cpp ) -src_use_openmp(equilibrium_stagger.cpp) -src_use_openmp(kinetic_stagger.cpp) - set( reactmicp_saturated_headers variablesfwd.hpp diffusive_upscaling_stagger.hpp equilibrium_stagger.hpp init_variables.hpp kinetic_stagger.hpp transport_program.hpp transport_stagger.hpp variables.hpp ) add_to_reactmicp_srcs_list(reactmicp_saturated_srcs) add_to_reactmicp_headers_list(reactmicp_saturated_headers) INSTALL(FILES ${reactmicp_saturated_headers} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/saturated_react ) diff --git a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp index b3bf6c9..dad8cd0 100644 --- a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp +++ b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp @@ -1,191 +1,190 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "equilibrium_stagger.hpp" #include "variables.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp_common/log.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" -#ifdef SPECMICP_USE_OPENMP -#include -#endif // SPECMICP_USE_OPENMP +#include "specmicp_common/config.h" + namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { using TrueConstPtr = SaturatedVariables * const; inline TrueConstPtr cast_to_var(solver::VariablesBase * const var) { return static_cast(var); } // EquilibriumStagger:: //! \brief Initialize the stagger at the beginning of an iteration void EquilibriumStagger::initialize_timestep( scalar_t dt, VariablesBase * const var) { m_dt = dt; TrueConstPtr true_var = cast_to_var(var); // Initialize velocity using values from previous timestep for (index_t node=0; nodenb_nodes(); ++node) { if (true_var->is_fixed_composition(node)) continue; scalar_t alpha = 1.0; for (index_t component: true_var->get_database()->range_aqueous_component()) { alpha = std::max(alpha, 0.9*dt*true_var->solid_concentration(node, component, true_var->chemistry_rate()) /(true_var->solid_concentration(node, component, true_var->displacement())) ); } auto solid_velocity = true_var->velocity().segment( true_var->offset_node(node)+true_var->offset_solid_concentration(), true_var->nb_component()); auto solid_chemistry_rate = true_var->chemistry_rate().segment( true_var->offset_node(node)+true_var->offset_solid_concentration(), true_var->nb_component()); solid_velocity = 1/alpha * solid_chemistry_rate; } } //! \brief Solve the equation for the timestep solver::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBase * const var) { TrueConstPtr true_var = cast_to_var(var); int failed_chemistry = 0; -#ifdef SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel default(none) shared(failed_chemistry) // node true_var being const is shared by default, can't be mentionned again { #pragma omp for schedule(dynamic, 5) for (index_t node=0; nodenb_nodes(); ++node) { // only solve if necessary if (true_var->is_fixed_composition(node) or failed_chemistry > 0) continue; const auto retcode = solve_one_node(node, true_var); if (retcode > 0) { ++failed_chemistry; } } } #else { for (index_t node=0; nodenb_nodes(); ++node) { if (true_var->is_fixed_composition(node)) continue; const auto retcode = solve_one_node(node, true_var); if (retcode > 0) { ++failed_chemistry; break; } } } -#endif // SPECMICP_USE_OPENMP +#endif // SPECMICP_HAVE_OPENMP if (failed_chemistry > 0) return solver::StaggerReturnCode::UnknownError; return solver::StaggerReturnCode::ResidualMinimized; } //! int EquilibriumStagger::solve_one_node( index_t node, SaturatedVariables * const true_var ) { AdimensionalSystemConstraints constraints(get_constraints(node)); constraints.total_concentrations = true_var->total_concentrations(node); AdimensionalSystemSolver adim_solver(true_var->get_database(), constraints, true_var->equilibrium_solution(node), m_options); Vector variables(true_var->equilibrium_solution(node).main_variables); micpsolver::MiCPPerformance perf = adim_solver.solve(variables); micpsolver::MiCPSolverReturnCode retcode = perf.return_code; if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) { ERROR << "Failed to solve chemistry problem at node " << node << ", return code = " << static_cast(retcode) << ", residual = " << perf.current_residual; ERROR << "Total concentration : \n" << constraints.total_concentrations; return 1; } true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables); AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node), true_var->get_database(), m_options.units_set); for (index_t component=0; componentnb_component(); ++component) { const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component); true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq; const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor()); const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt; true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq; const scalar_t c_sol = extractor.total_immobile_concentration(component); true_var->solid_concentration(node, component, true_var->displacement()) = c_sol; const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor()); const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt; true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol; true_var->solid_concentration(node, component, true_var->chemistry_rate()) = vel_sol; } return 0; } } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/CMakeLists.txt b/src/reactmicp/systems/unsaturated/CMakeLists.txt index 7078222..e791858 100644 --- a/src/reactmicp/systems/unsaturated/CMakeLists.txt +++ b/src/reactmicp/systems/unsaturated/CMakeLists.txt @@ -1,46 +1,43 @@ set( reactmicp_unsaturated_srcs aqueous_equation.cpp aqueous_pressure_equation.cpp boundary_conditions.cpp equilibrium_stagger.cpp pressure_equation.cpp saturation_equation.cpp saturation_pressure_equation.cpp transport_stagger.cpp variables.cpp variables_interface.cpp ) -src_use_openmp(equilibrium_stagger.cpp) -src_use_openmp(transport_stagger.cpp) - set( reactmicp_unsaturated_headers fv_1dof_equation.hpp fv_1dof_equation.inl init_variables.hpp simulation_data.hpp transport_constraints.hpp types_fwd.hpp upscaling_stagger_factory.hpp variables_box.hpp variables_sub.hpp aqueous_equation.hpp aqueous_pressure_equation.hpp boundary_conditions.hpp equilibrium_stagger.hpp pressure_equation.hpp saturation_equation.hpp saturation_pressure_equation.hpp transport_stagger.hpp variables.hpp variables_interface.hpp ) add_to_reactmicp_srcs_list(reactmicp_unsaturated_srcs) add_to_reactmicp_headers_list(reactmicp_unsaturated_headers) INSTALL(FILES ${reactmicp_unsaturated_headers} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/unsaturated ) diff --git a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp index fed3e8b..8381137 100644 --- a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp @@ -1,529 +1,530 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "equilibrium_stagger.hpp" #include "variables.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp_common/compat.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp_common/log.hpp" +#include "specmicp_common/config.h" namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { // ============================================ // // Declaration of implementation details // // ============================================ EquilibriumOptionsVector::EquilibriumOptionsVector(): utils::NameCachedVector() {} EquilibriumOptionsVector::EquilibriumOptionsVector( size_type size, std::string name, AdimensionalSystemSolverOptions& value ): utils::NameCachedVector( size, name, std::forward(value)) {} EquilibriumOptionsVector::EquilibriumOptionsVector( size_type size, std::string name ): utils::NameCachedVector(size, name) { static_assert(std::is_default_constructible::value, "AdimensionalSystemOptions must be default constructible."); } struct EquilibriumStagger::EquilibriumStaggerImpl { scalar_t m_dt {-1}; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_database; std::shared_ptr m_bcs; std::shared_ptr m_options; EquilibriumStaggerImpl( std::shared_ptr variables, std::shared_ptr boundary_conditions, std::shared_ptr options ): m_mesh(variables->get_mesh()), m_database(variables->get_database()), m_bcs(boundary_conditions), m_options(options) {} scalar_t dt() {return m_dt;} int compute_one_node(index_t node, UnsaturatedVariables * const vars); void compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ); void analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ); void initialize_timestep_one_node(index_t node, UnsaturatedVariables* vars); void initialize(UnsaturatedVariables* vars); void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars); StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); units::UnitsSet& get_units() { return m_options->get("default").units_set; } }; using TrueConstPtr = UnsaturatedVariables * const; inline TrueConstPtr cast_to_var(solver::VariablesBase * const var) { return static_cast(var); } // ============================================ // // Implementation // // ============================================ EquilibriumStagger::EquilibriumStagger( std::shared_ptr variables, std::shared_ptr boundary_conditions, std::shared_ptr options ): m_impl(utils::make_pimpl( variables, boundary_conditions, options)) { } EquilibriumStagger::~EquilibriumStagger() = default; //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize(true_vars); } //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize_timestep(dt, true_vars); } //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables EquilibriumStagger::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); return m_impl->restart_timestep(true_vars); } int EquilibriumStagger::EquilibriumStaggerImpl::compute_one_node( index_t node, UnsaturatedVariables * const vars ) { AdimensionalSystemConstraints constraints = m_bcs->get_constraint(node); compute_total_concentrations(node, constraints.total_concentrations, vars); // set partial pressure model { user_model_saturation_f callable = vars->get_vapor_pressure_model(); water_partial_pressure_f pv_model = [callable, node](scalar_t sat){ return callable(node, sat); }; constraints.set_water_partial_pressure_model(pv_model); } // We use the current value of the saturation to avoid converging // to the previous solution // This is necessary because the update to the saturation is very small const AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node); Vector x; AdimensionalSystemSolver the_solver; if (likely(previous_solution.is_valid)) // should always be true { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, previous_solution, m_options->operator [](node) ); x = previous_solution.main_variables; // we use the new value of the saturation to avoid a change // to small to be detected by the speciaion solver x(0) = porosity*saturation; // x(0) = volume fraction water } else { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, m_options->operator [](node) ); the_solver.initialise_variables(x, porosity*saturation, -3); } // micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { ERROR << "Failed to solve equilibrium at node " << node << ". \n" << "Return code : " << (int) perf.return_code; return -1; } AdimensionalSystemSolution solution = the_solver.get_raw_solution(x); analyse_solution(node, solution, vars); vars->set_adim_solution(node, solution); return 0; } void EquilibriumStagger::EquilibriumStaggerImpl::compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ) { total_concentrations.resize(m_database->nb_component()); const scalar_t saturation = vars->get_liquid_saturation()(node); const scalar_t porosity = vars->get_porosity()(node); const scalar_t rt = vars->get_rt(); // water { const scalar_t ctilde_w = vars->get_water_aqueous_concentration()(node); const scalar_t cbar_w = vars->get_solid_concentration(0)(node); scalar_t c_w = saturation*porosity*ctilde_w + cbar_w; if (vars->component_has_gas(0)) { const scalar_t pv_w = vars->get_pressure_main_variables(0)(node); c_w += (1.0-saturation)*porosity*pv_w/rt; } total_concentrations(0) = c_w; } total_concentrations(1) = 0.0; // aqueous components for (index_t component: m_database->range_aqueous_component()) { const scalar_t ctilde_i = vars->get_aqueous_concentration(component)(node); const scalar_t cbar_i = vars->get_solid_concentration(component)(node); scalar_t c_i = saturation*porosity*ctilde_i + cbar_i; if (vars->component_has_gas(component)) { const scalar_t pv_i = vars->get_pressure_main_variables(component)(node); c_i += (1.0-saturation)*porosity*pv_i/rt; } total_concentrations(component) = c_i; } } void EquilibriumStagger::EquilibriumStaggerImpl::analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ) { AdimensionalSystemSolutionExtractor extr(solution, m_database, get_units()); const scalar_t saturation = extr.saturation_water(); const scalar_t sat_g = 1.0 - saturation; const scalar_t rho_l = extr.density_water(); const scalar_t porosity = extr.porosity(); const scalar_t rt = vars->get_rt(); // porosity SecondaryTransientVariable& porosity_vars = vars->get_porosity(); const scalar_t porosity_0 = porosity_vars.predictor(node); const scalar_t porosity_vel = (porosity - porosity_0)/m_dt; porosity_vars.variable(node) = porosity; porosity_vars.velocity(node) = porosity_vel; // water // saturation MainVariable& satvars = vars->get_liquid_saturation(); const scalar_t saturation_0 = satvars.predictor(node); const scalar_t sat_vel = (saturation - saturation_0)/m_dt; satvars.variable(node) = saturation; satvars.velocity(node) = sat_vel; { // total aqueous concentration SecondaryTransientVariable& cwtilde_vars = vars->get_water_aqueous_concentration(); const scalar_t cwtilde = rho_l*extr.total_aqueous_concentration(0); const scalar_t cwtilde_0 = cwtilde_vars.predictor(node); const scalar_t cwtilde_vel = (cwtilde - cwtilde_0)/m_dt; cwtilde_vars.variable(node) = cwtilde; cwtilde_vars.velocity(node) = cwtilde_vel; // total solid concentration MainVariable& solid_vars = vars->get_solid_concentration(0); const scalar_t cwbar = extr.total_solid_concentration(0); const scalar_t cwbar_0 = solid_vars.predictor(node); const scalar_t cwbar_vel = (cwbar - cwbar_0)/m_dt; solid_vars.variable(node) = cwbar; solid_vars.velocity(node) = cwbar_vel; solid_vars.chemistry_rate(node) = - cwbar_vel; // vapor pressure if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); const scalar_t pv = vars->get_vapor_pressure_model()(node, saturation); const scalar_t pv_0 = pres_vars.predictor(node); const scalar_t pv_vel = (pv - pv_0)/m_dt; pres_vars.variable(node) = pv; pres_vars.velocity(node) = pv_vel; const scalar_t transient = ( (pv*porosity*sat_g) - (pv_0*porosity_0*(1.0-saturation_0)) )/ (rt*m_dt); const scalar_t pv_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pv_chem_rate; } } // aqueous components for (index_t component: m_database->range_aqueous_component()) { // total aqueous concentration MainVariable& aqconc = vars->get_aqueous_concentration(component); const scalar_t ctildei = rho_l*extr.total_aqueous_concentration(component); const scalar_t ctildei_0 = aqconc.predictor(node); const scalar_t ctildei_vel = (ctildei - ctildei_0) / m_dt; aqconc.variable(node) = ctildei; aqconc.velocity(node) = ctildei_vel; // total solid concentration MainVariable& solconc = vars->get_solid_concentration(component); const scalar_t cbari = extr.total_solid_concentration(component); const scalar_t cbari_0 = solconc.predictor(node); const scalar_t cbari_vel = (cbari - cbari_0) / m_dt; solconc.variable(node) = cbari; solconc.velocity(node) = cbari_vel; solconc.chemistry_rate(node) = - cbari_vel; // partial pressure if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); index_t id_g = vars->get_id_gas(component); const scalar_t pi = extr.fugacity_gas(id_g)*vars->get_total_pressure(); const scalar_t pi_0 = pres_vars.predictor(node); const scalar_t pi_vel = (pi - pi_0) / m_dt; pres_vars.variable(node) = pi; pres_vars.velocity(node) = pi_vel; const scalar_t transient = ( (pi*porosity*sat_g) - (pi_0*porosity_0*(1.0-saturation_0)) ) / (rt * m_dt); const scalar_t pi_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pi_chem_rate; } } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep_one_node( index_t node, UnsaturatedVariables * const vars ) { } EquilibriumStagger::StaggerReturnCode EquilibriumStagger::EquilibriumStaggerImpl::restart_timestep( UnsaturatedVariables* vars ) { int sum_retcode = 0; -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum_retcode) #endif for (index_t node=0; nodenb_nodes(); ++node) { if (not m_bcs->is_normal_node(node)) continue; sum_retcode += compute_one_node(node, vars); } if (sum_retcode != 0) { return StaggerReturnCode::UnknownError; } return StaggerReturnCode::ResidualMinimized; } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables* vars ) { m_dt = dt; // initialize //water { MainVariable& cbar_w = vars->get_solid_concentration(0); cbar_w.predictor = cbar_w.variable; cbar_w.velocity.setZero(); cbar_w.chemistry_rate.setZero(); SecondaryTransientVariable& ctilde_w = vars->get_water_aqueous_concentration(); ctilde_w.predictor = ctilde_w.variable; ctilde_w.velocity.setZero(); if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // aqueous component for (index_t component: m_database->range_aqueous_component()) { MainVariable& cbar_i = vars->get_solid_concentration(component); cbar_i.predictor = cbar_i.variable; cbar_i.velocity.setZero(); cbar_i.chemistry_rate.setZero(); if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // porosity { SecondaryTransientVariable& porosity = vars->get_porosity(); porosity.predictor = porosity.variable; porosity.velocity.setZero(); } for (auto node: m_mesh->range_nodes()) { initialize_timestep_one_node(node, vars); } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize( UnsaturatedVariables * const vars ) { // do nothing by default } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/transport_stagger.cpp b/src/reactmicp/systems/unsaturated/transport_stagger.cpp index 852df2c..9039d66 100644 --- a/src/reactmicp/systems/unsaturated/transport_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -1,1095 +1,1096 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "transport_stagger.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "variables.hpp" #include "variables_box.hpp" #include "transport_constraints.hpp" #include "saturation_equation.hpp" #include "saturation_pressure_equation.hpp" #include "pressure_equation.hpp" #include "aqueous_equation.hpp" #include "aqueous_pressure_equation.hpp" #include "specmicp_database/database_holder.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "specmicp_common/log.hpp" +#include "specmicp_common/config.h" #include namespace specmicp { namespace dfpmsolver { extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { //! \brief Type of the equation //! \internal enum class EquationType { Saturation, LiquidAqueous, Pressure }; //! \brief Format type to a string (for message error) //! \internal static std::string to_string(EquationType eq_type); //! \brief Format DFPM solver retcode to a string static std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode); // forward declaration of wrapper classes // They are defined at the bootom of this file class TaskBase; class SaturationTask; using VariablesBase = solver::VariablesBase; // =================================== // // // // Declaration of the implementation // // // // =================================== // //! \brief Implementation class for the Unsaturated transport stagger //! \internal //! //! This class does all the work //! It was designed to be OpenMP compatible class UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl { public: // Constructor // ----------- UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ); // Main functions // -------------- //! \brief Return the residual scalar_t get_residual(UnsaturatedVariables * const vars); //! \brief Return the first residual (start of timestep) scalar_t get_residual_0(UnsaturatedVariables * const vars); //! \brief Return the update (norm of velocity vector) scalar_t get_update(UnsaturatedVariables * const vars); //! \brief Initialize timestep void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars); //! \brief Restart the timestep solver::StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); // Timestep management // ------------------- //! \brief Save the timestep void save_dt(scalar_t dt) {m_dt = dt;} //! \brief Return the timestep scalar_t get_dt() {return m_dt;} // Task index // ---------- index_t& id_aqueous_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } index_t id_aqueous_task(index_t component) const { return m_aq_equation_task[static_cast(component)]; } index_t& id_gas_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } index_t id_gas_tasj(index_t component) const { return m_aq_equation_task[static_cast(component)]; } //! \brief Return the options of the saturation equation dfpmsolver::ParabolicDriverOptions* get_saturation_options(); //! \brief Return the options of aqueous equation of component dfpmsolver::ParabolicDriverOptions* get_aqueous_options(index_t component); //! \brief Return the options of gas equation of component dfpmsolver::ParabolicDriverOptions* get_gas_options(index_t component); private: scalar_t m_norm_0 {1.0}; // timestep management scalar_t m_dt {-1.0}; //! \brief The saturation equations std::unique_ptr m_saturation_equation; // Equation and solver std::vector> m_equation_list; std::vector m_aq_equation_task; std::vector m_gas_equation_task; std::vector m_merged_gas; }; // Main functions // =============== // call of the implementation UnsaturatedTransportStagger::UnsaturatedTransportStagger( std::shared_ptr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_impl(utils::make_pimpl( variables, boundary_conditions, merge_saturation_pressure, merge_aqueous_pressure )) { } UnsaturatedTransportStagger::~UnsaturatedTransportStagger() = default; static UnsaturatedVariables * const get_var(VariablesBase * const var) { return static_cast(var); } void UnsaturatedTransportStagger::initialize_timestep( scalar_t dt, VariablesBase* var ) { m_impl->initialize_timestep(dt, get_var(var)); } solver::StaggerReturnCode UnsaturatedTransportStagger::restart_timestep(VariablesBase * const var) { return m_impl->restart_timestep(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual(VariablesBase * const var) { return m_impl->get_residual(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual_0(VariablesBase * const var) { return m_impl->get_residual_0(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_update(VariablesBase * const var) { return m_impl->get_update(get_var(var)); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_saturation_options() { return m_impl->get_saturation_options(); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_aqueous_options(index_t component) { return m_impl->get_aqueous_options(component); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_gas_options(index_t component) { return m_impl->get_gas_options(component); } // =============================== // // =============================== // // // // Implementation details // // ---------------------- // // // // =============================== // // =============================== // // 2 main sections // - Solver wrappers : wrapper around 1 couple equation/solver // - UnsaturatedTransportStaggerImpl : call the wrappers // =============================== // // // // Solver wrappers // // // // =============================== // // This wrappers are used to abstract the residual computation // and the methods to solve every equations //! \brief Base class for a task //! //! A task is how we solve governing equations, //! and obtain residuals and update //! //! \internal class SPECMICP_DLL_LOCAL TaskBase { public: TaskBase(index_t component, EquationType eq_type): m_type(eq_type), m_component(component) {} virtual ~TaskBase() {} //! \brief Compute the squared residuals virtual scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared residuals at the beginning of the timestep virtual scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared update of the variables virtual scalar_t compute_squared_update( UnsaturatedVariables * const vars ) = 0; //! \brief Initialize the timestep virtual void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) = 0; //! \brief Solve the governing equation virtual dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) = 0; //! \brief Return the component index (in the db) index_t component() {return m_component;} //! \brief The equation type EquationType equation_type() {return m_type;} private: EquationType m_type; index_t m_component; }; //! \brief Base class for a equation solver //! //! \internal template class SPECMICP_DLL_LOCAL EquationTask: public TaskBase { public: using EqT = typename traits::EqT; using SolverT = typename traits::SolverT; EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, const TransportConstraints& constraints ): TaskBase(component, traits::equation_type), m_equation(the_mesh, variables, constraints), m_solver(m_equation) {} EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, const TransportConstraints& constraints, scalar_t scaling ): TaskBase(component, traits::equation_type), m_equation(the_mesh, variables, constraints), m_solver(m_equation) { m_equation.set_scaling(scaling); } Derived* derived() {return static_cast(this);} // This function must be implemented by subclass MainVariable& get_var(UnsaturatedVariables * const vars) { return derived()->get_var(vars); } scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, true); return residuals.squaredNorm()/m_norm_square_0; } scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, false); m_norm_square_0 = residuals.squaredNorm(); if (m_norm_square_0 < 1e-8) m_norm_square_0 = 1.0; return 1.0; } scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { const MainVariable& main = get_var(vars); return main.velocity.squaredNorm(); } void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) override { m_dt = dt; MainVariable& main = get_var(vars); main.predictor = main.variable; main.transport_fluxes.setZero(); main.velocity.setZero(); m_solver.initialize_timestep(dt, get_var(vars).variable); } dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) override { MainVariable& main = get_var(vars); m_solver.velocity() = main.velocity; auto retcode = m_solver.restart_timestep(main.variable); if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet) { main.velocity = m_solver.velocity(); } return retcode; } dfpmsolver::ParabolicDriverOptions& get_options() { return m_solver.get_options(); } const dfpmsolver::ParabolicDriverOptions& get_options() const { return m_solver.get_options(); } scalar_t get_dt() {return m_dt;} protected: scalar_t m_norm_square_0 {-1}; scalar_t m_dt {-1}; EqT m_equation; SolverT m_solver; }; //! \brief Traits struct for the SaturationTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationTaskTraits { using EqT = SaturationEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationTask: public EquationTask { using base = EquationTask; public: SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints ): EquationTask( component, the_mesh, variables, constraints) {} SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints, scalar_t scaling ): EquationTask( component, the_mesh, variables, constraints, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the SaturationPressureTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationPressureTaskTraits { using EqT = SaturationPressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationPressureVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationPressureTask: public EquationTask { using base = EquationTask; public: SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints ): EquationTask( component, the_mesh, variables, constraints) {} SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints, scalar_t scaling ): EquationTask( component, the_mesh, variables, constraints, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidAqueousTaskTraits { using EqT = AqueousTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidAqueousComponentVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidAqueousTask: public EquationTask { using base = EquationTask; public: LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints ): EquationTask( component, the_mesh, variables, constraints) {} LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints, scalar_t scaling ): EquationTask( component, the_mesh, variables, constraints, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidGasAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidGasAqueousTaskTraits { using EqT = AqueousGasTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidGasAqueousVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidGasAqueousTask: public EquationTask { using base = EquationTask; public: LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints ): base(component, the_mesh, variables, constraints) {} LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints, scalar_t scaling ): base(component, the_mesh, variables, constraints, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the Pressure Task traits //! //! \internal struct SPECMICP_DLL_LOCAL PressureTaskTraits { using EqT = PressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = PressureVariableBox; static constexpr EquationType equation_type {EquationType::Pressure}; }; //! \brief Wrapper for the pressure equation solver //! //! \internal class SPECMICP_DLL_LOCAL PressureTask: public EquationTask { using base = EquationTask; public: PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints ): EquationTask( component, the_mesh, variables, constraints) {} PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, const TransportConstraints& constraints, scalar_t scaling ): EquationTask( component, the_mesh, variables, constraints, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_pressure_main_variables(component()); } }; // ================================== // // // // UnsaturatedTransportStaggerImpl // // // // ================================== // // constructor // =========== UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_merged_gas(variables->get_database()->nb_component(), false) { database::RawDatabasePtr raw_db = variables->get_database(); mesh::Mesh1DPtr the_mesh = variables->get_mesh(); m_aq_equation_task = std::vector( raw_db->nb_component(), no_equation); m_gas_equation_task = std::vector( raw_db->nb_component(), no_equation); TransportConstraints constraints; constraints.set_fixed_nodes(boundary_conditions->get_fixed_nodes()); constraints.set_gas_nodes(boundary_conditions->get_gas_nodes()); dfpmsolver::ParabolicDriverOptions* opts = nullptr; // Saturation equations // ==================== // // There is 2 choices for the saturation equation // Either SaturationPressureEquation or SaturationPressure if (merge_saturation_pressure and variables->component_has_gas(0)) { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_pressure_variables(), constraints, variables->get_aqueous_scaling(0) ); m_merged_gas[0] = true; opts = &(static_cast( m_saturation_equation.get())->get_options()); } else { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_variables(), constraints, variables->get_aqueous_scaling(0) ); opts = &(static_cast( m_saturation_equation.get())->get_options()); } opts->residuals_tolerance = 1e-8; opts->absolute_tolerance = 1e-12; opts->step_tolerance = 1e-12; opts->threshold_stationary_point = 1e-12; opts->maximum_iterations = 500; opts->sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts->quasi_newton = 3; opts->maximum_step_length = 0.1; opts->max_iterations_at_max_length = 5000; const index_t size = raw_db->nb_aqueous_components() + variables->nb_gas(); m_equation_list.reserve(size); // Liquid aqueous diffusion-advection // ================================== for (index_t id: raw_db->range_aqueous_component()) { if (merge_aqueous_pressure and variables->component_has_gas(id)) { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_gas_aqueous_variables(id), constraints, variables->get_aqueous_scaling(id)) ); m_merged_gas[id] = true; } else { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_aqueous_component_variables(id), constraints, variables->get_aqueous_scaling(id)) ); } id_aqueous_task(id) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_aqueous_options(id); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } // Water partial pressure // ====================== // // Depending on the saturation equation chosen we may or may not include // this equations if (variables->component_has_gas(0) and (not m_merged_gas[0])) { m_equation_list.emplace_back( make_unique(0, the_mesh, variables->get_pressure_variables(0), constraints, variables->get_gaseous_scaling(0) )); id_gas_task(0) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_gas_options(0); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } // Gaseous diffusion for (index_t id: raw_db->range_aqueous_component()) { if (variables->component_has_gas(id) and (not m_merged_gas[id])) { m_equation_list.emplace_back( make_unique(id, the_mesh, variables->get_pressure_variables(id), constraints, variables->get_gaseous_scaling(id) )); id_gas_task(id) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_gas_options(id); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } } } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_saturation_options() { dfpmsolver::ParabolicDriverOptions* opts = nullptr; if (m_merged_gas[0]) { opts = &static_cast( m_saturation_equation.get())->get_options(); } else { opts = &static_cast( m_saturation_equation.get())->get_options(); } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_aqueous_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_aqueous_task(component); if (id != no_equation) { if (m_merged_gas[component]) { opts = &static_cast( m_equation_list[id].get())->get_options(); } else { opts = &static_cast( m_equation_list[id].get())->get_options(); } } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_gas_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_gas_task(component); if (id != no_equation) { opts = &static_cast( m_equation_list[id].get())->get_options(); } return opts; } // Residuals // ========= scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual( UnsaturatedVariables * const vars ) { auto sum = m_saturation_equation->compute_squared_residual(vars); -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_residual(vars); } auto norm = std::sqrt(sum); //std::cout << "; " << norm << std::endl; return norm; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual_0( UnsaturatedVariables * const vars ) { return m_norm_0; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_update( UnsaturatedVariables * const vars ) { scalar_t sum = m_saturation_equation->compute_squared_update(vars); -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_update(vars); } return std::sqrt(sum); } // Solving the equations // ===================== void UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) { save_dt(dt); vars->get_advection_flux().set_zero(); m_saturation_equation->initialize_timestep(dt, vars); scalar_t sum = m_saturation_equation->compute_squared_residual_0(vars); -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqinitialize_timestep(dt, vars); sum += task->compute_squared_residual_0(vars); } m_norm_0 = 1.0;//std::sqrt(sum); } solver::StaggerReturnCode UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep( UnsaturatedVariables * const vars ) { dfpmsolver::ParabolicDriverReturnCode retcode; bool flag_fail = false; bool flag_error_minimized = false; // Saturation equation retcode = m_saturation_equation->restart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Failed to solve saturation equation, return code :" << to_string(retcode); flag_fail = true; } if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized) flag_error_minimized = true; // Other equation if (not flag_fail) { vars->set_relative_variables(); -#if SPECMICP_USE_OPENMP +#ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(or: flag_fail) \ reduction(and: flag_error_minimized) #endif for (std::size_t ideq=0; ideqrestart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Equation of type '" << to_string(task->equation_type()) << "' for component " << task->component() << " has failed with return code : " << to_string(retcode); flag_fail = true; } flag_error_minimized = flag_error_minimized and (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized); } } // Return code solver::StaggerReturnCode return_code = solver::StaggerReturnCode::ResidualMinimized; if (flag_error_minimized) { return_code = solver::StaggerReturnCode::ErrorMinimized; } if (flag_fail) return_code = solver::StaggerReturnCode::UnknownError; return return_code; } // ================================ // // // // Helper functions // // // // ================================ // static std::string to_string(EquationType eq_type) { std::string str; switch (eq_type) { case EquationType::LiquidAqueous: str = "Liquid advection-diffusion"; break; case EquationType::Pressure: str = "Gaseous diffusion"; break; case EquationType::Saturation: str = "Saturation equation"; break; default: str = "Unknown equation"; break; } return str; } std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode) { using RetCode = dfpmsolver::ParabolicDriverReturnCode; std::string str; switch (retcode) { case RetCode::ResidualMinimized: str = "ResidualMinimized"; break; case RetCode::ErrorMinimized: str = "ErrorMinimized"; break; case RetCode::NotConvergedYet: str = "NotConvergedYet"; break; case RetCode::StationaryPoint: str = "StationaryPoint"; break; case RetCode::ErrorLinearSystem: str = "ErrorLinearSystem"; break; case RetCode::MaxStepTakenTooManyTimes: str = "MaxStepTakenTooManyTimes"; break; case RetCode::MaxIterations: str = "MaxIterations"; break; case RetCode::LinesearchFailed: str = "LinesearchFailed"; default: str = "Unknown return code"; break; } return str; } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/specmicp_common/CMakeLists.txt b/src/specmicp_common/CMakeLists.txt index 3e236d7..175a385 100644 --- a/src/specmicp_common/CMakeLists.txt +++ b/src/specmicp_common/CMakeLists.txt @@ -1,133 +1,138 @@ # Main # ==== # var to store the files set(specmicp_common_srcs "" CACHE INTERNAL "specmicp_common files" FORCE) set(specmicp_common_headers "" CACHE INTERNAL "specmicp_common headers" FORCE ) # macro to add to vars macro(add_to_main_srcs_list LIST_NAME) set(tmp "") foreach (src ${${LIST_NAME}}) list(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${src}) endforeach(src) set(specmicp_common_srcs "${specmicp_common_srcs};${tmp}" CACHE INTERNAL "specmicp common files" FORCE) endmacro(add_to_main_srcs_list) macro(add_to_main_headers_list LIST_NAME) set( tmp "") foreach(header ${${LIST_NAME}}) LIST(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${header}) endforeach(header) set(specmicp_common_headers "${specmicp_common_headers};${tmp}" CACHE INTERNAL "headers for specmicp common" FORCE) endmacro(add_to_main_headers_list) # set( specmicp_common_main_srcs dateandtime.cpp filesystem.cpp log.cpp moving_average.cpp string_algorithms.cpp timer.cpp ) add_to_main_srcs_list( specmicp_common_main_srcs ) set( specmicp_common_main_headers cached_vector.hpp compat.hpp dateandtime.hpp filesystem.hpp log.hpp macros.hpp moving_average.hpp options_handler.hpp perfs_handler.hpp pimpl_ptr.hpp range_iterator.hpp range_iterator.inl scope_guard.hpp string_algorithms.hpp timer.hpp types.hpp ) add_to_main_headers_list( specmicp_common_main_headers ) INSTALL(FILES ${headers} DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp_common ) add_subdirectory(cli) add_subdirectory(eigen) add_subdirectory(io) add_subdirectory(micpsolver) add_subdirectory(odeint) add_subdirectory(physics) add_subdirectory(plugins) add_subdirectory(sparse_solvers) +set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS config.h.in) +add_custom_target(config_h SOURCES config.h.in) +configure_file(config.h.in + ${CMAKE_CURRENT_BINARY_DIR}/config.h + ) # build library # ============= set_pgo_flag(${specmicp_common_srcs}) add_library(objspecmicp_common OBJECT ${specmicp_common_srcs} ${specmicp_common_headers} ) set_property(TARGET objspecmicp_common PROPERTY POSITION_INDEPENDENT_CODE 1) add_library(specmicp_common SHARED $) if (UNIX) set(specmicp_common_link_libraries dl) else() message(FATAL_ERROR "Plugin system only for POSIX at this time !") endif() if (HDF5_FOUND) set(specmicp_common_link_libraries "${HDF5_LIBRARIES};${specmicp_common_link_libraries}") endif() target_link_libraries(specmicp_common ${specmicp_common_link_libraries}) install(TARGETS specmicp_common LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(specmicp_common_static STATIC $ ) install(TARGETS specmicp_common_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(specmicp_common_static EXCLUDE_FROM_ALL STATIC $ ) endif() target_link_libraries(specmicp_common_static ${specmicp_common_link_libraries}) set_target_properties(specmicp_common_static PROPERTIES OUTPUT_NAME specmicp_common) diff --git a/src/specmicp_common/config.h.in b/src/specmicp_common/config.h.in new file mode 100644 index 0000000..9501bae --- /dev/null +++ b/src/specmicp_common/config.h.in @@ -0,0 +1,47 @@ +/* ============================================================================= + + Copyright (c) 2014 - 2016 + F. Georget Princeton University + All rights reserved. + + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * + +============================================================================= */ + +#ifndef SPECMICP_COMMON_CONFIG_H +#define SPECMICP_COMMON_CONFIG_H + +//! \file config.h.in +//! \brief Compilation configuration header + +// defined if has the secure_getenv +#cmakedefine SPECMICP_HAVE_SECURE_GETENV + +// define if openmp is set +#cmakedefine SPECMICP_HAVE_OPENMP + +#endif // SPECMICP_COMMON_CONFIG_H diff --git a/src/specmicp_common/filesystem.cpp b/src/specmicp_common/filesystem.cpp index 8b435cf..713c47d 100644 --- a/src/specmicp_common/filesystem.cpp +++ b/src/specmicp_common/filesystem.cpp @@ -1,211 +1,213 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "filesystem.hpp" #include #include #include #include #include #include #include +#include "specmicp_common/config.h" + #include "log.hpp" static std::string test_name; namespace specmicp { namespace utils { bool is_directory(const std::string& path) { struct stat info; if (stat(path.c_str(), &info) == -1) { return false; } if (S_ISDIR(info.st_mode)) { return true; } return false; } bool is_file(const std::string& path) { struct stat info; if (stat(path.c_str(), &info) == -1) { return false; } if (S_ISREG(info.st_mode)) { return true; } return false; } std::string get_current_directory() { char buf[PATH_MAX]; char* test = getcwd(buf, PATH_MAX); if (test == NULL) { ERROR << "Unable to obtain current working directory"; throw std::runtime_error("Something is wrong," "unable to obtain the current directory."); } return std::string(buf); } std::string complete_path( const std::string& dir, const std::string& file ) { std::string complete = dir; if (dir.back() != '/') { complete += '/'; } complete += file; return complete; } std::string relative_to_absolute( const std::string& rel_path, std::string& error ) { std::string abs_path = ""; char buf[PATH_MAX]; // call posix function char* res = realpath(rel_path.c_str(), buf); if (res == NULL) { // parse error if (errno == ENOENT) { error = "No such file '" + rel_path +"'."; } else if (errno == EACCES) { error = "Read permission denied while searching for '" + rel_path +"'."; } else if (errno == EIO) { error = "I/O error while searching for '" + rel_path +"'."; } else { error = "Error while accessing '" + rel_path + "'."; } } else { // no error, copy buffer abs_path = buf; } return abs_path; } int name_filter(const struct dirent64* entry) { auto res = strncmp(entry->d_name, test_name.c_str(), 256); if (res == 0) { return 1; } return 0; } std::string find_path( std::string filename, const std::vector &directories ) { // check if filename is a path const auto has_sep = filename.find('/'); if (has_sep != std::string::npos) { // already a path => we convert it to absolute std::string error = ""; auto filepath = relative_to_absolute(filename, error); if (filepath == "") { ERROR << error; return ""; // empty string is signal for error } return filepath; } // if not a path we try to find it std::string complete_path_str = ""; test_name = filename.c_str(); for (auto dir: directories) { bool found = false; struct dirent64** entry_list; auto count = scandir64(dir.c_str(), &entry_list, name_filter, alphasort64); if (count < 0) { ERROR << "Problem while scanning directory : " << dir << "."; throw std::runtime_error("Problem while scanning directory " + dir + "."); } if (count == 0) { continue; } if (count > 1) { WARNING << "More that one match for file '" << filename << "in : " << dir << "."; } for (auto ind=0; indd_name); } free(entry); // need to free everything } free(entry_list); if (found) break; } return complete_path_str; } std::string get_env( const std::string& env_var ) { -#ifdef _GNU_SOURCE +#ifdef SPECMICP_HAVE_SECURE_GETENV char* env = secure_getenv(env_var.c_str()); #else char* env = getenv(env_var.c_str()); #endif if (env == NULL) { return ""; } else return env; } } // end namespace utils } // end namespace specmicp