diff --git a/CMakeLists.txt b/CMakeLists.txt index 130ccbb..48a9060 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,352 +1,354 @@ project(specmicp) cmake_minimum_required(VERSION 2.8.8) ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package find_package(CxxTest) # Necessary if you want the tests # catch for unit testing include(ExternalProject) find_package(Git REQUIRED) ExternalProject_Add( catch PREFIX ${CMAKE_BINARY_DIR}/catch GIT_REPOSITORY https://github.com/philsquared/Catch.git TIMEOUT 10 UPDATE_COMMAND ${GIT_EXECUTABLE} pull CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" LOG_DOWNLOAD ON ) ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### set(PROJECT_TEST_DIR ${PROJECT_SOURCE_DIR}/tests) ExternalProject_Get_Property(catch source_dir) set(CATCH_INCLUDE_DIR ${source_dir}/include CACHE INTERNAL "Path to include folder for Catch") set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${EIGEN3_INCLUDE_DIR}) include_directories(${PROJECT_SOURCE_DIR}) include_directories(${CATCH_INCLUDE_DIR}) include_directories(${PROJECT_TEST_DIR}) ################ Modules ########################## # -------- MiCPSolver ---------------- set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(specmisolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolverold.hpp ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # ------- SpecMiCP --------------- set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/simulation_box.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ${SPECMICP_DIR}/speciation_solver.cpp ${SPECMICP_DIR}/simple_composition.cpp ) 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 ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # --------- ODE int ------------------ set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ) # --------- DFPMSolver ----------------- set(DFPMSOLVER_DIR ${PROJECT_SOURCE_DIR}/dfpmsolver) add_custom_target(dfpmsolver_incl SOURCES ${DFPMSOLVER_DIR}/driver.hpp ${DFPMSOLVER_DIR}/driver.inl ${DFPMSOLVER_DIR}/driver_structs.hpp ${DFPMSOLVER_DIR}/dfpm_program.hpp # ${DFPMSOLVER_DIR}/sparse_solver.hpp ${DFPMSOLVER_DIR}/sparse_solver_structs.hpp # Parabolic solver ${DFPMSOLVER_DIR}/parabolic_driver.hpp ${DFPMSOLVER_DIR}/parabolic_driver.inl ${DFPMSOLVER_DIR}/parabolic_program.hpp ${DFPMSOLVER_DIR}/parabolic_structs.hpp ) # --------- ReactMiCP ----------------- set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_DIR}/common.hpp # ${REACTMICP_DIR}/micpsolvers/parabolicprogram.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver_structs.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.inl # ${REACTMICP_DIR}/meshes/mesh1d.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.inl ${REACTMICP_DIR}/systems/saturated_diffusion/options.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.hpp + ${REACTMICP_DIR}/systems/saturated_diffusion/boundary_conditions.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_DIR}/systems/diffusion/diffusion.cpp ${REACTMICP_DIR}/systems/diffusion/diffusion_numbering.cpp ${REACTMICP_DIR}/systems/diffusion/diffusion_secondary.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/equation_numbering.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/speciation_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/variables.cpp ) add_library(reactmicp ${REACTMICP_LIBFILES}) # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp ) # ------- Physics ---------------- set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # ------- Data ------------------- set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # ------ Documentation ----------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) endif(DOXYGEN_FOUND) # --------- Python interface ---------- add_subdirectory( cython ) # --------- Test ---------------------- enable_testing(true) # CXX Test # ======== if(CXXTEST_FOUND) set(CXXTEST_USE_PYTHON TRUE) include_directories(${CXXTEST_INCLUDE_DIR}) enable_testing() # MiCP Solver # ------------ CXXTEST_ADD_TEST(test_cond_number test_cond_number.cpp ${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp) CXXTEST_ADD_TEST(test_ncp_funtion test_ncp_function.cpp ${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp) CXXTEST_ADD_TEST(test_micpsolver test_micpsolver.cpp ${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp) # Database # -------- CXXTEST_ADD_TEST(test_data_reader test_data_reader.cpp ${PROJECT_TEST_DIR}/database/test_data_reader.hpp) target_link_libraries(test_data_reader specmicp_database jsoncpp) CXXTEST_ADD_TEST(test_data_selector test_data_selector.cpp ${PROJECT_TEST_DIR}/database/test_data_selector.hpp) target_link_libraries(test_data_selector specmicp_database jsoncpp) # Odeint # ------ CXXTEST_ADD_TEST(test_butchertableau test_butchertableau.cpp ${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp) CXXTEST_ADD_TEST(test_embeddedrungekutta test_embeddedrungekutta.cpp ${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp) endif() # Test 'Perso' # ============ add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp) target_link_libraries(thermocarbo specmicp specmicp_database jsoncpp) add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp) target_link_libraries(carboalu specmicp specmicp_database jsoncpp) add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp) target_link_libraries(alkaliactivated specmicp specmicp_database jsoncpp) add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp) target_link_libraries(rate_nicoleau specmicp specmicp_database jsoncpp) add_executable(kinetics ${PROJECT_TEST_DIR}/specmicp/kinetics.cpp) target_link_libraries(kinetics specmicp specmicp_database jsoncpp) set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv ) add_custom_target(data_nicoleau SOURCES ${DATA_TEST}) file(INSTALL ${DATA_TEST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/) # Catch test #=========== add_custom_target(catch_header SOURCES ${PROJECT_TEST_DIR}/test_database.hpp ) add_executable(test_reactmicp_diffusion ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/test_reactmicp_diffusion.cpp ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/diffusion.cpp ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/numbering.cpp ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/secondary.cpp ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/solving.cpp ${PROJECT_TEST_DIR}/reactmip/systems/diffusion/utils.cpp ) target_link_libraries(test_reactmicp_diffusion reactmicp specmicp specmicp_database jsoncpp) add_executable(test_reactmicp_saturated_diffusion ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/test_reactmicp_saturated_diffusion.cpp ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/variables.cpp ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/solver.cpp ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/transport_program.cpp + ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/reactive_transport.cpp ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/utils.cpp ) target_link_libraries(test_reactmicp_saturated_diffusion reactmicp specmicp specmicp_database jsoncpp) diff --git a/src/reactmicp/systems/saturated_diffusion/boundary_conditions.hpp b/src/reactmicp/systems/saturated_diffusion/boundary_conditions.hpp new file mode 100644 index 0000000..f09e7ef --- /dev/null +++ b/src/reactmicp/systems/saturated_diffusion/boundary_conditions.hpp @@ -0,0 +1,44 @@ +#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP +#define SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP + +#include +#include "specmicp/equilibrium_data.hpp" + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace siasaturated { + +//! \brief The different types of Boundary conditions accepted +enum class BoundaryConditionType +{ + Free, //!< No BC + FixedComposition //!< Fixed composition for this node +}; + +//! \brief Store the boundary conditions to be parsed by the reactive transport solver +struct SIABoundaryConditions +{ + //! \brief list of the different initial state used + //! + //! The first one should be the one used everywhere + std::vector list_initial_states; + //! \brief The BC type for each node + std::vector bs_types; + //! \brief The initial state to use for each node + std::vector initial_states; + + SIABoundaryConditions(int nb_nodes): + bs_types(nb_nodes, BoundaryConditionType::Free), + initial_states(nb_nodes, 0) + {} +}; + + +} // end namespace siasaturated +} // end namespace systems +} // end namespace reactmicp +} // end namespace specmicp + + +#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/equation_numbering.hpp b/src/reactmicp/systems/saturated_diffusion/equation_numbering.hpp deleted file mode 100644 index b48c0a3..0000000 --- a/src/reactmicp/systems/saturated_diffusion/equation_numbering.hpp +++ /dev/null @@ -1,4 +0,0 @@ -#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP -#define SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP - -#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_EQUATIONNUMBERING_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp index 3051823..b9c715b 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp @@ -1,7 +1,25 @@ #include "reactive_transport_solver.hpp" namespace specmicp { namespace reactmicp { + +void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs) +{ + assert(bcs.bs_types.size() == m_mesh->nnodes()); + for (int node: m_mesh->range_nodes()) + { + switch (bcs.bs_types[node]) { + case BoundaryConditionType::FixedComposition: + m_bcs(node) = 1; + break; + default: + m_bcs(node) = 0; + break; + } + m_variables.update_composition(node, bcs.list_initial_states[bcs.initial_states[node]]); + } +} + } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp index 990eb7f..3c65bb8 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp @@ -1,43 +1,72 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP #include "options.hpp" #include "variables.hpp" +#include "transport_program.hpp" +#include "dfpmsolver/parabolic_driver.hpp" +#include "boundary_conditions.hpp" namespace specmicp { namespace reactmicp { using namespace systems::siasaturated; //! \brief SIA solver for saturated reactive transport problem //! class SIASaturatedReactiveTransportSolver { public: + SIASaturatedReactiveTransportSolver( + std::shared_ptr the_mesh, + std::shared_ptr the_database, + std::shared_ptr transport_parameters + ): + m_mesh(the_mesh), + m_variables(the_database, the_mesh, transport_parameters), + m_transport_program(the_mesh, the_database, transport_parameters), + m_transport_solver(m_transport_program), + m_bcs(the_mesh->nnodes()) + {} + + void apply_boundary_conditions(SIABoundaryConditions bcs); + + SIASaturatedVariables& get_variables() {return m_variables;} + //! \brief Solve a timestep void solve_timestep(double dt); void restart_transport_iterations(); void restart_speciation_iterations(); void solve_speciation_node(int node); void check_convergence(); //! \brief Return a read-write reference to the options SIASaturatedReactiveTransportSolverOptions& get_options() {return m_options;} //! \brief Return a read-only reference to the options const SIASaturatedReactiveTransportSolverOptions& get_options() const {return m_options;} private: void initialize(); SIASaturatedReactiveTransportSolverOptions m_options; + std::shared_ptr m_mesh; + SIASaturatedVariables m_variables; + + // transport program + SaturatedDiffusionProgram m_transport_program; + dfpmsolver::ParabolicDriver m_transport_solver; + + // + Eigen::VectorXi m_bcs; + }; } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/variables.cpp b/src/reactmicp/systems/saturated_diffusion/variables.cpp index 4193b00..6ad779c 100644 --- a/src/reactmicp/systems/saturated_diffusion/variables.cpp +++ b/src/reactmicp/systems/saturated_diffusion/variables.cpp @@ -1,59 +1,59 @@ #include "variables.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { void SIASaturatedVariables::update_composition(int node, EquilibriumState& composition) { // speciation m_speciation_variables.col(node) = composition.main_variables(); m_secondary_concentration.col(node) = composition.molality_secondary(); m_secondary_variables.col(node) = composition.activity_coefficient(); m_ionic_strength(node) = composition.ionic_strength(); // transport for (int component: m_data->range_aqueous_component()) { total_mobile_concentration(node, component) = composition.total_aqueous_concentration_component(component); } - m_speciation_variables(node, 0) = mass_water(node); + m_speciation_variables(0, node) = mass_water(node); } double SIASaturatedVariables::immobile_total_amount(int node, int component) { double sum = 0; for (int mineral: m_data->range_mineral()) { if (m_data->nu_mineral(mineral, component) != 0) { sum += m_data->nu_mineral(mineral, component)*mineral_amount(node, mineral); } } return sum; } void SIASaturatedVariables::nodal_update_total_amount(int node, Eigen::VectorXd& total_amount) { total_amount.resize(m_data->nb_component); total_amount(0) = mass_water(node)/m_data->molar_mass_basis_si(0); for (int component: m_data->range_aqueous_component()) { total_amount(component) = nodal_component_update_total_amount(node, component); } } double SIASaturatedVariables::nodal_component_update_total_amount(int node, int component) { return mass_water(node)*total_mobile_concentration(node,component)+immobile_total_amount(node, component); } double SIASaturatedVariables::mass_water(int node) { return m_param->density_water()*m_param->porosity(node)*m_mesh->cell_volume(node); } } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/variables.hpp b/src/reactmicp/systems/saturated_diffusion/variables.hpp index 35dc510..4a7fd6b 100644 --- a/src/reactmicp/systems/saturated_diffusion/variables.hpp +++ b/src/reactmicp/systems/saturated_diffusion/variables.hpp @@ -1,141 +1,140 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP #include "reactmicp/systems/secondary_variables/secondary.hpp" #include "specmicp/equilibrium_data.hpp" #include "reactmicp/meshes/mesh1d.hpp" #include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { //! Variables for the SIASaturatedReactiveTransportSolver class SIASaturatedVariables: public SecondaryVariables { public: SIASaturatedVariables( - int nb_nodes, std::shared_ptr thedatabase, std::shared_ptr themesh, std::shared_ptr theparam ): - SecondaryVariables(nb_nodes, thedatabase), + SecondaryVariables(themesh->nnodes(), thedatabase), m_data(thedatabase), m_mesh(themesh), m_param(theparam), m_total_mobile_concentration( - Eigen::VectorXd((thedatabase->nb_component-1)* nb_nodes)), + Eigen::VectorXd((thedatabase->nb_component-1)* themesh->nnodes())), m_speciation_variables( - Eigen::MatrixXd(thedatabase->nb_component+thedatabase->nb_mineral, nb_nodes)) + Eigen::MatrixXd(thedatabase->nb_component+thedatabase->nb_mineral, themesh->nnodes())) {} int get_dof_total_concentration(int node, int component) const { assert(component >= 1 and component < m_data->nb_component); - assert(node > 0 and node < m_speciation_variables.cols()); + assert(node >= 0 and node < m_speciation_variables.cols()); return (component+offset_total_component())+(m_data->nb_component+offset_total_component())*node; } //! Return the total mobile concentration of 'component' at 'node' double total_mobile_concentration(int node, int component) const { return m_total_mobile_concentration.coeff(get_dof_total_concentration(node, component)); } //! Return a reference to the total mobile concentration of 'component' at 'node' double& total_mobile_concentration(int node, int component) { return m_total_mobile_concentration.coeffRef(get_dof_total_concentration(node, component)); } Eigen::VectorXd& total_mobile_concentrations() { return m_total_mobile_concentration; } Eigen::Block total_mobile_concentrations(int node) { return m_total_mobile_concentration.segment(node*(m_data->nb_component+offset_total_component()), (m_data->nb_component+offset_total_component())); } //! Return the concentration of 'component' at 'node' double component_concentration(int node, int component) const { assert(component >= 1 and component < m_data->nb_component); assert(node >= 0 and node < m_speciation_variables.cols()); return pow10(m_speciation_variables.coeff(component+offset_component(), node)); } //! Return the concentration of 'component' at 'node' double log_component_concentration(int node, int component) const { assert(component >= 1 and component < m_data->nb_component); assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeff(component+offset_component(), node); } //! Return a reference to the concentration of 'component' at 'node' double& log_component_concentration(int node, int component) { assert(component >= 1 and component < m_data->nb_component); assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeffRef(component+offset_component(), node); } //! Return the number of moles of 'mineral' at 'node' double mineral_amount(int node, int mineral) const { assert(mineral >= 0 and mineral < m_data->nb_mineral); assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeff(offset_mineral()+mineral, node); } //! Return a reference to the number of moles of 'mineral' at 'node' double& mineral_amount(int node, int mineral) { assert(mineral >= 0 and mineral < m_data->nb_mineral); assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeffRef(offset_mineral()+mineral, node); } //! Reference to the speciation variables Eigen::MatrixXd& speciation_variables() { return m_speciation_variables; } //! Eigen Expression to the speciation variables Eigen::MatrixXd::ColXpr speciation_variables(int node) { return m_speciation_variables.col(node); } //! Return the chemical composition of a node - needed to initialize the speciation solver EquilibriumState equilibrium_composition(int node) { return EquilibriumState(m_speciation_variables.col(node), m_secondary_concentration.col(node), m_secondary_variables.col(node), m_ionic_strength(node), m_data ); } //! Update the composition at 'node' using the EquilibriumState instance 'composition' void update_composition(int node, EquilibriumState& composition); //! Return the total moles number for 'component' at 'node' from the speciation double immobile_total_amount(int node, int component); //! Return the mass of water at 'node' double mass_water(int node); //! Compute the total moles number at 'node' using the transport mobile concentration void nodal_update_total_amount(int node, Eigen::VectorXd& total_amount); private: // where the component informations are stored int offset_total_component() const {return -1;} // where the component informations are stored int offset_component() const {return 0;} // where the mineral informations are stored int offset_mineral() const {return m_data->nb_component+offset_component();} // Used in the update FEM->speciation double nodal_component_update_total_amount(int node, int component); std::shared_ptr m_data; std::shared_ptr m_mesh; std::shared_ptr m_param; Eigen::VectorXd m_total_mobile_concentration; Eigen::MatrixXd m_speciation_variables; }; } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP diff --git a/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp b/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp new file mode 100644 index 0000000..973c236 --- /dev/null +++ b/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp @@ -0,0 +1,46 @@ +#include "catch.hpp" + +#include "utils.hpp" +#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp" +#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp" + +using namespace specmicp; +using namespace specmicp::reactmicp; +using namespace specmicp::reactmicp::systems; +using namespace specmicp::reactmicp::systems::siasaturated; + + +TEST_CASE("reactive transport solver", "[reactive transport, speciation, dfpm, solver]") { + + int nb_nodes = 5; + std::shared_ptr database = get_test_carbo_database(); + + std::shared_ptr themesh = std::make_shared(nb_nodes-1, 1e-3, 0.001); + auto parameters = std::make_shared(1e-8, 0.2); + + SIABoundaryConditions bcs(nb_nodes); + bcs.list_initial_states.push_back(sample_carbo_composition(database)); + bcs.list_initial_states.push_back(blank_composition(database)); + + bcs.bs_types[0] = BoundaryConditionType::FixedComposition; + bcs.initial_states[0] = 1; + + + SECTION("Initialization") { + SIASaturatedReactiveTransportSolver solver(themesh, database, parameters); + + solver.apply_boundary_conditions(bcs); + + SIASaturatedVariables& variables = solver.get_variables(); + + REQUIRE(std::abs(variables.mineral_amount(1, 1) + - bcs.list_initial_states[0].moles_mineral(1)) < 1e-10 ); + REQUIRE(std::abs(variables.component_concentration(1, 1) + - bcs.list_initial_states[0].molality_component(1)) < 1e-10 ); + REQUIRE(std::abs(variables.mineral_amount(0, 1) + - bcs.list_initial_states[1].moles_mineral(1)) < 1e-10 ); + REQUIRE(std::abs(variables.component_concentration(0, 1) + - bcs.list_initial_states[1].molality_component(1)) < 1e-10 ); + + } +} diff --git a/tests/reactmip/systems/saturated_diffusion/solver.cpp b/tests/reactmip/systems/saturated_diffusion/solver.cpp index 61798e7..a2db7b3 100644 --- a/tests/reactmip/systems/saturated_diffusion/solver.cpp +++ b/tests/reactmip/systems/saturated_diffusion/solver.cpp @@ -1,88 +1,88 @@ #include "catch.hpp" #include "utils.hpp" #include "dfpmsolver/parabolic_driver.hpp" #include "reactmicp/systems/saturated_diffusion/transport_program.hpp" #include "reactmicp/systems/saturated_diffusion/variables.hpp" #include "specmicp/reduced_system_solver.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("dfpm solver", "[dfpm, solver]") { SECTION("transport program") { std::shared_ptr database = get_test_carbo_database(); int nb_nodes = 5; EquilibriumState initial_state = sample_carbo_composition(database); std::shared_ptr themesh = std::make_shared(nb_nodes-1, 1e-3, 0.001); auto parameters = std::make_shared(1e-8, 0.2); - SIASaturatedVariables variables(nb_nodes, database, themesh, parameters); + SIASaturatedVariables variables(database, themesh, parameters); SaturatedDiffusionProgram the_program(themesh, database, parameters); Eigen::VectorXd& concentrations = variables.total_mobile_concentrations(); concentrations.setOnes(); concentrations.block(0, 0, the_program.get_ndf(), 1).setZero(); the_program.number_equations({0,}); dfpmsolver::ParabolicDriver solver(the_program); solver.get_velocity() = Eigen::VectorXd::Zero(concentrations.rows()); double dt = 100; dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(dt, concentrations); REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); for (int node=1; nodennodes(); ++node) { for (int component: database->range_aqueous_component()) { the_program.external_flux(node, component) = parameters->porosity(0)*(1.0 - concentrations(the_program.get_dof_component(node, component))) / (parameters->density_water()*dt*themesh->cell_volume(node)); } } retcode = solver.restart_timestep(concentrations); REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); } SECTION("micpsolver") { std::shared_ptr database = get_test_carbo_database(); database::Database database_handler(database); int nb_nodes = 5; std::shared_ptr themesh = std::make_shared(nb_nodes-1, 1e-3, 0.001); auto parameters = std::make_shared(1e-8, 0.2); - SIASaturatedVariables variables(nb_nodes, database, themesh, parameters); + SIASaturatedVariables variables(database, themesh, parameters); EquilibriumState initial_state = sample_carbo_composition(database); variables.update_composition(1, initial_state); Eigen::VectorXd totconc; variables.nodal_update_total_amount(1, totconc); totconc(database_handler.component_label_to_id("HCO3[-]")) += 0.05; totconc(database_handler.component_label_to_id("HO[-]")) -= 0.05; ReducedSystemSolverOptions options; options.solver_options.max_factorization_step = 1; ReducedSystemSolver solver(database, totconc, options); Eigen::VectorXd spec_var(variables.speciation_variables().col(1)); specmicp::micpsolver::MiCPPerformance perf = solver.solve(spec_var); REQUIRE((int) perf.return_code == (int) specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized); EquilibriumState final_compo = solver.get_solution(spec_var); variables.update_composition(1, final_compo); } } diff --git a/tests/reactmip/systems/saturated_diffusion/variables.cpp b/tests/reactmip/systems/saturated_diffusion/variables.cpp index 8998a8e..f504cc9 100644 --- a/tests/reactmip/systems/saturated_diffusion/variables.cpp +++ b/tests/reactmip/systems/saturated_diffusion/variables.cpp @@ -1,64 +1,64 @@ #include "catch.hpp" #include "utils.hpp" #include "reactmicp/systems/saturated_diffusion/variables.hpp" #include "reactmicp/meshes/mesh1d.hpp" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::siasaturated; TEST_CASE("SIA Variables", "[Variables, Secondary variables, SIA]") { std::shared_ptr database = get_test_carbo_database(); int nb_nodes = 5; std::shared_ptr themesh = std::make_shared(nb_nodes-1, 1e-3, 0.001); auto parameters = std::make_shared(1e-6, 0.2); SECTION("Initialization") { - SIASaturatedVariables variables(nb_nodes, database, themesh, parameters); + SIASaturatedVariables variables(database, themesh, parameters); REQUIRE_NOTHROW( variables.component_concentration(0, 2) ); REQUIRE_NOTHROW( variables.component_concentration(2, 2) ); REQUIRE_NOTHROW( variables.component_concentration(4, 4) ); REQUIRE_NOTHROW( variables.component_concentration(3, 1) ); REQUIRE_NOTHROW( variables.log_component_concentration(3, 3) ); REQUIRE_NOTHROW( variables.log_gamma_component(4, 2) ); } SECTION("Ionic strength computation") { - SIASaturatedVariables variables(nb_nodes, database, themesh, parameters); + SIASaturatedVariables variables(database, themesh, parameters); database::Database database_handler(database); variables.log_component_concentration(1, database_handler.component_label_to_id("Ca[2+]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("HO[-]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("HCO3[-]")) = -3; variables.log_component_concentration(1, database_handler.component_label_to_id("SiO(OH)3[-]")) = -3; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3)*0.5 ); variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO3[2-]")) = 1e-2; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 ); variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO2")) = 1e-1; variables.nodal_ionic_strength(1); REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 ); } SECTION("Update") { - SIASaturatedVariables variables(nb_nodes, database, themesh, parameters); + SIASaturatedVariables variables(database, themesh, parameters); EquilibriumState composition = sample_carbo_composition(database); variables.update_composition(1, composition); REQUIRE(variables.component_concentration(1, 1) == composition.molality_component(1)); REQUIRE(variables.total_mobile_concentration(1, 1) == composition.total_aqueous_concentration_component(1)); } }