diff --git a/CMakeLists.txt b/CMakeLists.txt index c44112e..7e55728 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,487 +1,492 @@ #################################################### # # # SpecMiCP : CMakeLists.txt # # # #################################################### project(specmicp) cmake_minimum_required(VERSION 2.8.8) # For an explanation of the options see the INSTALL file option(SPECMICP_USE_OPENMP "Thread parallelisation of specmicp instances if available" ON) # External Package ######################################################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") # Boost # ----- find_package(Boost REQUIRED) include_directories(${Boost_INCLUDE_DIR}) # 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 if(EXISTS "${EIGEN3_INCLUDE_DIR}/unsupported/Eigen/src/IterativeSolvers/GMRES.h") add_definitions(-DEIGEN_UNSUPPORTED_FOUND) INCLUDE_DIRECTORIES("${EIGEN3_INCLUDE_DIR}/unsupported/") endif() # Compilation flags ######################################################################### # OpenMP # ------ include(CheckCXXCompilerFlag) check_cxx_compiler_flag("-fopenmp" HAS_OPENMP) # c++11 # ----- check_cxx_compiler_flag("-std=c++11" HAS_CXX11) if(NOT HAS_CXX11) message(FATAL_ERROR "A c++11 compatible compiler is necessary") endif() if(UNIX) # Use ld.gold if it is available and isn't disabled explicitly # Ref : https://bugs.webkit.org/show_bug.cgi?id=137953 option(USE_LD_GOLD "Use GNU gold linker" ON) if (USE_LD_GOLD) execute_process(COMMAND ${CMAKE_C_COMPILER} -fuse-ld=gold -Wl,--version ERROR_QUIET OUTPUT_VARIABLE LD_VERSION) if ("${LD_VERSION}" MATCHES "GNU gold") set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fuse-ld=gold") else () message(WARNING "GNU gold linker isn't available, using the default system linker.") endif () endif () # Generic options SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") if (HAS_OPENMP) SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") endif() SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -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() # Activate this option to use finite difference jacobian # add_definitions(-DSPECMICP_DEBUG_EQUATION_FD_JACOBIAN) # Directories ######################################################################### # 3rd party # ========= # JsonCPP # ------- set(JSONCPP_DIR ${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) add_custom_target(3rdparty_header SOURCES ${JSONCPP_DIR}/json/json.h ${JSONCPP_DIR}/json/json-forwards.h ) include_directories(${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) # Source # ======= set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # Includes # ========= include_directories(${PROJECT_SOURCE_DIR}) add_custom_target(common SOURCES ${PROJECT_SOURCE_DIR}/common.hpp ${PROJECT_SOURCE_DIR}/database.hpp ) # Modules ######################################################################### # MiCPSolver # ============== # Mixed complementarity Solver # include only module set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(micpsolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolverold.hpp ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # SpecMiCP # ============ # The speciation solver set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) set(KINETICS_DIR ${SPECMICP_DIR}/kinetics) set(ADIMENSIONAL_DIR ${SPECMICP_DIR}/adimensional) set(ADIMKINETICS_DIR ${SPECMICP_DIR}/adimensional_kinetics) set(PROBLEM_SOLVER_DIR ${SPECMICP_DIR}/problem_solver) # headers only module add_custom_target(specmic_inc SOURCES ${SPECMICP_DIR}/simulation_box.hpp ${SPECMICP_DIR}/reduced_system_structs.hpp ${SPECMICP_DIR}/reduced_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_numbering.hpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm_structs.hpp ${PROBLEM_SOLVER_DIR}/formulation.hpp ${KINETICS_DIR}/kinetic_model.hpp ${KINETICS_DIR}/kinetic_variables.hpp ${KINETICS_DIR}/kinetic_system_solver_structs.hpp ${ADIMKINETICS_DIR}/kinetic_model.hpp ${ADIMKINETICS_DIR}/kinetic_variables.hpp ${ADIMKINETICS_DIR}/kinetic_system_solver_structs.hpp ) # Source for the library set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ${ADIMENSIONAL_DIR}/adimensional_system.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solver.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solution_extractor.cpp ${ADIMENSIONAL_DIR}/adimensional_system_pcfm.cpp ${ADIMENSIONAL_DIR}/equilibrium_curve.cpp ${PROBLEM_SOLVER_DIR}/dissolver.cpp # kinetic ${KINETICS_DIR}/kinetic_system.cpp ${KINETICS_DIR}/kinetic_system_solver.cpp ${ADIMKINETICS_DIR}/kinetic_system.cpp ${ADIMKINETICS_DIR}/kinetic_system_solver.cpp ${ADIMKINETICS_DIR}/kinetic_system_euler_solver.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # Database # ============ set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/data_container.cpp ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ${JSONCPP_DIR}/jsoncpp.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # ODEint # ========== # ODE integration, headers-only set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ${ODEINT_DIR}/runge_kutta_step_structs.hpp ) # DFPMSolver # ============== ## Parabolic solver for finite element problem, header only set(DFPMSOLVER_DIR ${PROJECT_SOURCE_DIR}/dfpmsolver) add_custom_target(dfpmsolver_incl SOURCES ${DFPMSOLVER_DIR}/driver.hpp ${DFPMSOLVER_DIR}/driver.inl ${DFPMSOLVER_DIR}/driver_structs.hpp ${DFPMSOLVER_DIR}/dfpm_program.hpp # # Parabolic solver ${DFPMSOLVER_DIR}/parabolic_driver.hpp ${DFPMSOLVER_DIR}/parabolic_driver.inl ${DFPMSOLVER_DIR}/parabolic_program.hpp ${DFPMSOLVER_DIR}/parabolic_structs.hpp ) # DFPM # ======== # Dynaflow+- , Finite element solver, mostly headers set (DFPM_DIR ${PROJECT_SOURCE_DIR}/dfpm) add_custom_target(dfpm_incl SOURCES ${DFPM_DIR}/types.hpp ${DFPM_DIR}/1dtransport/diffusion_parameters.hpp ${DFPM_DIR}/mesh.hpp ${DFPM_DIR}/meshes/mesh1dfwd.hpp ${DFPM_DIR}/meshes/mesh1d.hpp ${DFPM_DIR}/meshes/uniform_mesh1d.hpp ${DFPM_DIR}/meshes/generic_mesh1d.hpp ${DFPM_DIR}/meshes/axisymmetric_uniform_mesh1d.hpp ${DFPM_DIR}/meshes/axisymmetric_mesh1d.hpp ) set(DFPMLIB ${DFPM_DIR}/1dtransport/diffusion.cpp ) add_library(dfpm ${DFPMLIB}) # ReactMiCP # ============= # Reactive transport solver set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp) set(REACTMICP_SOLVER_DIR ${REACTMICP_DIR}/solver/) set(REACTMICP_SYSTEMS_DIR ${REACTMICP_DIR}/systems/) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_DIR}/common.hpp # ${REACTMICP_DIR}/micpsolvers/parabolicprogram.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver_structs.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver.inl ## remove or fix # ${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.inl ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver_structs.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters_base.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/boundary_conditions.hpp # Flexible reactive tranport solver # --------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_variables.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/variables.cpp ${REACTMICP_DIR}/equilibrium_curve/chemistry.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_extractor.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_coupler.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ${REACTMICP_SOLVER_DIR}/timestepper.cpp # ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ) add_library(reactmicp ${REACTMICP_LIBFILES}) if(HAS_OPENMP) if(SPECMICP_USE_OPENMP) SET_TARGET_PROPERTIES(reactmicp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP") endif() endif() # Utils # ========= # Stuff that may be used by the other modules set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) set(IO_DIR ${UTILS_DIR}/io) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp # sparse solvers # -------------- ${UTILS_DIR}/sparse_solvers/sparse_solver.hpp ${UTILS_DIR}/sparse_solvers/sparse_solver_base.hpp ${UTILS_DIR}/sparse_solvers/sparse_solver_structs.hpp ${UTILS_DIR}/sparse_solvers/sparse_qr.hpp ${UTILS_DIR}/sparse_solvers/sparse_lu.hpp ${UTILS_DIR}/sparse_solvers/sparse_bicgstab.hpp ${UTILS_DIR}/sparse_solvers/sparse_gmres.hpp ${UTILS_DIR}/options_handler.hpp ${UTILS_DIR}/perfs_handler.hpp ${UTILS_DIR}/moving_average.hpp ${UTILS_DIR}/timer.hpp # input/output + + ${IO_DIR}/units.hpp + ${IO_DIR}/meshes.hpp ${IO_DIR}/reactive_transport.hpp + ${IO_DIR}/saturated_react.hpp + ) # Physics # =========== # Stuff that may be used by the other modules, but is related to physics set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # Database ######################################################################### # Not code, just necessary data set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ${DATA_DIR}/momas_benchmark.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # Documentation ######################################################################### # "common" documentation # ----------------------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # Doxygen documentation # --------------------- # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) # Configure doxygen configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) # Citations for the documentations 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/) # Target to build the documentations 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 ) endif(DOXYGEN_FOUND) # The following modules have their own CMakeLists.txt # Python interface ######################################################################### add_subdirectory( cython ) # Tests ######################################################################### add_subdirectory( tests ) # Examples ######################################################################### add_subdirectory( examples ) diff --git a/src/utils/io/meshes.hpp b/src/utils/io/meshes.hpp new file mode 100644 index 0000000..4791c4a --- /dev/null +++ b/src/utils/io/meshes.hpp @@ -0,0 +1,34 @@ +#ifndef SPECMICP_IO_MESHES_HPP +#define SPECMICP_IO_MESHES_HPP + +#include +#include +#include "dfpm/meshes/mesh1d.hpp" +#include "physics/units.hpp" +#include "utils/io/units.hpp" + +namespace specmicp { +namespace io { + +inline void print_mesh(std::ostream* output, mesh::Mesh1DPtr the_mesh, units::LengthUnit lenght_unit) +{ + (*output) << "# dfpm mesh \n"; + (*output) << "# units : " << length_unit_to_string(lenght_unit) << "\n"; + (*output) << "Node\tDistance\tCell_volume" << std::endl; + for (index_t node: the_mesh->range_nodes()) + { + (*output) << node << "\t" << the_mesh->get_position(node) << "\t" << the_mesh->get_volume_cell(node) << "\n"; + } +} + +inline void print_mesh(std::string filepath, mesh::Mesh1DPtr the_mesh, units::LengthUnit lenght_unit) +{ + std::ofstream outfile(filepath); + print_mesh(&outfile, the_mesh, lenght_unit); + outfile.close(); +} + +} // end namespace io +} // end namespace specmicp + +#endif // SPECMICP_IO_MESHES_HPP diff --git a/src/utils/io/reactive_transport.hpp b/src/utils/io/reactive_transport.hpp index fdf6b53..236856e 100644 --- a/src/utils/io/reactive_transport.hpp +++ b/src/utils/io/reactive_transport.hpp @@ -1,139 +1,161 @@ #ifndef SPECMICP_IO_REACTMICP_HPP #define SPECMICP_IO_REACTMICP_HPP #include #include "reactmicp/solver/reactive_transport_solver_structs.hpp" +#include "utils/timer.hpp" + namespace specmicp { namespace io { //! \brief Print the header inline void print_reactmicp_header( std::ostream* output ); //! \brief Print the performance information of the reactive transport solver inline void print_reactmicp_performance_short( std::ostream* output, const reactmicp::solver::ReactiveTransportPerformance& perf ); //! \brief Print the timing information of the reactive transport solver inline void print_reactmicp_timer( std::ostream* output, const reactmicp::solver::ReactiveTransportTimer& timer ); //! \brief Return a literal version of a ReactiveTransportReturnCode inline std::string reactmicp_return_code_to_string( reactmicp::solver::ReactiveTransportReturnCode retcode); inline void print_reactmicp_header( std::ostream* output ) { (*output) << "// ================================== //\n" << "// //\n" << "// ReactMiCP //\n" << "// //\n" << "// ================================== //\n\n" << "A reactive transport solver based on SpecMiCP.\n" << "------------------------------------------- \n" << std::endl; } inline void print_reactmicp_performance_short( std::ostream* output, const reactmicp::solver::ReactiveTransportPerformance& perf ) { (*output) << "------------------------------------------- \n" << " Performance \n " << "------------------------------------------- \n\t" << " - Return code : " << reactmicp_return_code_to_string(perf.return_code) << "\n\t" << " - Residuals : " << perf.residuals << "\n\t" << " - Number of iterations : " << perf.nb_iterations << "\n" << "------------------------------------------- " << std::endl; } inline void print_reactmicp_timer( std::ostream* output, const reactmicp::solver::ReactiveTransportTimer& timer ) { (*output) << "------------------------------------------- \n" - << "Time spend in each stagger \n" + << "Time spent in each stagger \n" << "------------------------------------------- \n\t" << "- Transport : " << timer.transport_time << " s\n\t" << "- Chemistry : " << timer.chemistry_time << " s\n\t" << "- Upscaling : " << timer.upscaling_time << " s\n" << "-----------------------------------------" << std::endl; } +inline void print_reactmicp_end( + std::ostream* output, + const Timer& total_timer, + const reactmicp::solver::ReactiveTransportTimer& timer + ) +{ + scalar_t elapsed_s = total_timer.elapsed_time(); + index_t hours = static_cast(elapsed_s) / 3600; + index_t elapsed_minutes = elapsed_s - hours*3600; + index_t minute = elapsed_minutes / 60; + index_t seconds = elapsed_minutes - 60*minute; + + (*output) << " ====================================================== \n"; + (*output) << "computation finished at " + << total_timer.get_ctime_stop(); + (*output) << " Duration of the computation : " << total_timer.elapsed_time() << " s" + << "( " << hours << "h " << minute << "min " << seconds << "s )\n"; + print_reactmicp_timer(output, timer); +} + inline void print_reactmicp_performance_long_header( std::ostream* output ) { (*output) << "Id\t T\t dt\t Return_code\t Iterations \t Residuals\t Time\t Transport_time\t Chemistry_time" << std::endl; } inline void print_reactmicp_performance_long( std::ostream* output, index_t cnt, scalar_t total, const reactmicp::solver::ReactiveTransportPerformance& perfs ) { (*output) << cnt << "\t " << total << "\t " << perfs.timestep << "\t " << (int) perfs.return_code << "\t " << perfs.nb_iterations << "\t " << perfs.residuals << "\t " << perfs.total_time << "\t" << perfs.transport_time << "\t" << perfs.chemistry_time << std::endl; } inline std::string reactmicp_return_code_to_string( reactmicp::solver::ReactiveTransportReturnCode retcode) { using RetCode = reactmicp::solver::ReactiveTransportReturnCode; switch (retcode) { case RetCode::ResidualMinimized: return "ResidualMinimized"; case RetCode::ErrorMinimized: return "ErrorMinimized"; default: switch (retcode) { case RetCode::StationaryPoint: return "StationaryPoint"; case RetCode::MaximumIterationsReached: return "MaximumIterationsReached"; case RetCode::GoodEnough: return "Good enough..."; case RetCode::TransportBypass: return "Transport Bypass"; case RetCode::TransportFailure: return "TransportFailure"; case RetCode::ChemistryFailure: return "ChemistryFailure"; case RetCode::UpscalingFailure: return "UpscalingFailure"; default: return "Unknow error code !"; } } } } // end namespace io } // end namespace specmicp #endif // SPECMICP_IO_REACTMICP_HPP diff --git a/src/utils/io/saturated_react.hpp b/src/utils/io/saturated_react.hpp new file mode 100644 index 0000000..267e923 --- /dev/null +++ b/src/utils/io/saturated_react.hpp @@ -0,0 +1,168 @@ +#ifndef SPECMICP_IO_SATURATEDREACT_HPP +#define SPECMICP_IO_SATURATEDREACT_HPP + +#include +#include + +#include "utils/io/units.hpp" + +//#include "reactmicp/systems/saturated_react/variablesfwd.hpp" +#include "reactmicp/systems/saturated_react/variables.hpp" +#include +#include + +using namespace specmicp::reactmicp::systems::satdiff; + +namespace specmicp { +namespace io { + +static constexpr char field_sep = '\t'; +static constexpr char comment = '#'; + +void print_csv_header(std::ostream* ofile) +{ + std::time_t time = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now()); + (*ofile) << comment << " - ReactMiCP - " << std::ctime(&time); +} + +void print_components_total_aqueous_concentration( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, + std::ostream* ofile + ) +{ + print_csv_header(ofile); + (*ofile) << comment << "units : mol/" << io::volume_unit_to_string(units_set.length) << std::endl; + (*ofile) << "Position"; + for (index_t component: the_database->range_aqueous_component()) + { + (*ofile) << field_sep << the_database->labels_basis[component]; + } + (*ofile) << std::endl; + + for (index_t node: the_mesh->range_nodes()) + { + (*ofile) << the_mesh->get_position(node); + + for (index_t component: the_database->range_aqueous_component()) + { + (*ofile) << field_sep << variables->aqueous_concentration(node, component); + } + (*ofile) << std::endl; + } +} + +void print_components_total_aqueous_concentration( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, + const std::string& filepath + ) +{ + std::ofstream ofile(filepath); + print_components_total_aqueous_concentration(the_database, variables, the_mesh, units_set, &ofile); + ofile.close(); +} + +void print_components_total_solid_concentration( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, + std::ostream* ofile + ) +{ + print_csv_header(ofile); + (*ofile) << comment << "units : mol/" << io::volume_unit_to_string(units_set.length) << std::endl; + (*ofile) << "Position"; + for (index_t component: the_database->range_aqueous_component()) + { + (*ofile) << field_sep << the_database->labels_basis[component]; + } + (*ofile) << std::endl; + + for (index_t node: the_mesh->range_nodes()) + { + (*ofile) << the_mesh->get_position(node); + + for (index_t component: the_database->range_aqueous_component()) + { + (*ofile) << field_sep << variables->solid_concentration(node, component); + } + (*ofile) << std::endl; + } +} + +void print_components_total_solid_concentration( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + const units::UnitsSet& units_set, + const std::string& filepath + ) +{ + std::ofstream ofile(filepath); + print_components_total_solid_concentration(the_database, variables, the_mesh, units_set, &ofile); + ofile.close(); +} + +void print_minerals_profile( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + std::ostream* ofile + ) +{ + print_csv_header(ofile); + (*ofile) << "Position"; + for (index_t mineral: the_database->range_mineral()) + { + (*ofile) << field_sep << the_database->labels_minerals[mineral]; + } + (*ofile) << field_sep << "Porosity"; + (*ofile) << field_sep << "pH"; + (*ofile) << std::endl; + + for (index_t node: the_mesh->range_nodes()) + { + (*ofile) << the_mesh->get_position(node); + AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet()); + + for (index_t mineral: the_database->range_mineral()) + { + (*ofile) << field_sep << extractor.volume_fraction_mineral(mineral); + } + (*ofile) << field_sep << variables->porosity(node); + (*ofile) << field_sep << extractor.pH(); + (*ofile) << std::endl; + + } + +} + +void print_minerals_profile( + RawDatabasePtr the_database, + SaturatedVariablesPtr variables, + mesh::Mesh1DPtr the_mesh, + const std::string& filepath + ) +{ + std::ofstream ofile(filepath); + print_minerals_profile( + the_database, + variables, + the_mesh, + &ofile + ); + ofile.close(); +} + + + +} // end namespace io +} // end namespace specmicp + +#endif // SPECMICP_IO_SATURATEDREACT_HPP diff --git a/src/utils/io/units.hpp b/src/utils/io/units.hpp new file mode 100644 index 0000000..d691b3d --- /dev/null +++ b/src/utils/io/units.hpp @@ -0,0 +1,34 @@ +#ifndef SPECMICP_IO_UNITS_HPP +#define SPECMICP_IO_UNITS_HPP + +#include +#include "physics/units.hpp" + +namespace specmicp { +namespace io { + +std::string length_unit_to_string(units::LengthUnit length_u) +{ + switch (length_u) { + case units::LengthUnit::meter: + return "m"; + break; + case units::LengthUnit::decimeter: + return "dm"; + break; + case units::LengthUnit::centimeter: + return "cm"; + break; + } +} + +std::string volume_unit_to_string(units::LengthUnit length_u) +{ + return length_unit_to_string(length_u) + "^3"; +} + + +} // end namespace io +} // end namespace specmicp + +#endif // SPECMICP_IO_UNITS_HPP