diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index ba990e2..5c990e4 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,153 +1,159 @@ # ReactMiCP # ============= set(REACTMICP_SOLVER_DIR solver) set(REACTMICP_SYSTEMS_DIR systems) set(REACTMICP_EQUILIBRIUMCURVE equilibrium_curve) add_custom_target(reactmicp_incl SOURCES ${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 + + ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_sub.hpp + ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_box.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ${REACTMICP_SOLVER_DIR}/timestepper.cpp ${REACTMICP_SOLVER_DIR}/runner.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/kinetic_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.cpp + + ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables.cpp + io/reactive_transport.cpp io/saturated_react.cpp io/configuration.cpp ) set_visibility_flag(${REACTMICP_LIBFILES}) set_pgo_flag(${REACTMICP_LIBFILES}) if(OPENMP_FOUND) if(SPECMICP_USE_OPENMP) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) endif() endif() add_library(objreactmicp OBJECT ${REACTMICP_LIBFILES}) set_property(TARGET objreactmicp PROPERTY POSITION_INDEPENDENT_CODE 1) add_library(reactmicp SHARED $) target_link_libraries(reactmicp dfpm) install(TARGETS reactmicp LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # include files # ------------- set(REACTMICP_STAGGERBASE_INCLUDE_LIST ${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 ) set(REACTMICP_SOLVER_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.hpp ${REACTMICP_SOLVER_DIR}/timestepper.hpp ${REACTMICP_SOLVER_DIR}/runner.hpp ) set(REACTMICP_SATURATEDREACT_INCLUDE_LIST ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.hpp ) set(REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp ) set(REACTMICP_IO_INCLUDE_LIST io/reactive_transport.hpp io/saturated_react.hpp io/configuration.hpp ) install(FILES ${REACTMICP_STAGGERBASE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver/staggers_base ) install(FILES ${REACTMICP_SOLVER_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver ) install(FILES ${REACTMICP_SATURATEDREACT_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/saturated_react ) install(FILES ${REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/equilibrium_curve ) install(FILES ${REACTMICP_IO_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/io ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(reactmicp_static STATIC $) install(TARGETS reactmicp_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(reactmicp_static EXCLUDE_FROM_ALL STATIC $) endif() set_target_properties(reactmicp_static PROPERTIES OUTPUT_NAME reactmicp) target_link_libraries(reactmicp_static dfpm_static) diff --git a/src/reactmicp/systems/unsaturated/variables.cpp b/src/reactmicp/systems/unsaturated/variables.cpp new file mode 100644 index 0000000..6f7698b --- /dev/null +++ b/src/reactmicp/systems/unsaturated/variables.cpp @@ -0,0 +1,248 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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 "variables.hpp" +#include "../../../dfpm/meshes/mesh1d.hpp" + +#include "variables_box.hpp" + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + +// =================== +// List main variables +// =================== + +ListMainVariable::ListMainVariable(index_t nb_vars, index_t nb_nodes) +{ + variables.reserve(nb_vars); + variables.emplace_back(make_unique(nb_nodes)); + variables.emplace_back(nullptr); + for (index_t i=2; i(nb_nodes)); + } +} + +ListMainVariable::ListMainVariable( + index_t nb_vars, + index_t nb_nodes, + std::vector to_init + ) +{ + variables.reserve(nb_vars); + for (const auto& is_present: to_init) { + if (is_present) { + variables.emplace_back(make_unique(nb_nodes)); + } else { + variables.emplace_back(nullptr); + } + } +} + +//! \brief Reset the variables +void ListMainVariable::reset() +{ + variables[0]->reset(); + for (vector_type::size_type i=2; ireset(); + } +} + + +void ListMainVariable::hard_reset(index_t component, index_t nb_nodes) +{ + specmicp_assert(component >= 0 and (unsigned) component < variables.size()); + if (nb_nodes == 0) + variables[component].reset(nullptr); + else + variables[component].reset(new MainVariable(nb_nodes)); +} + + + +// ============== +// Variables box +// ============== + +SaturationVariableBox::SaturationVariableBox(UnsaturatedVariables& vars): + liquid_saturation(vars.m_liquid_var.water()), + solid_concentration(vars.m_solid_var.water()), + vapor_pressure(vars.m_gas_var.water()), + aqueous_concentration(vars.m_water_aq_concentration), + porosity(vars.m_porosity), + liquid_permeability(vars.m_liquid_permeability), + liquid_diffusivity(vars.m_liquid_diffusivity), + relative_liquid_permeability(vars.m_relative_liquid_permeability), + relative_liquid_diffusivity(vars.m_relative_liquid_diffusivity), + capillary_pressure(vars.m_capillary_pressure), + advection_flux(vars.m_advection_flux), + capillary_pressure_f(vars.m_capillary_pressure_f), + relative_liquid_permeability_f(vars.m_relative_liquid_permeability_f), + relative_liquid_diffusivity_f(vars.m_relative_liquid_diffusivity_f) +{} + +LiquidAqueousComponentVariableBox::LiquidAqueousComponentVariableBox(UnsaturatedVariables& vars, index_t component): + aqueous_concentration(vars.m_liquid_var.aqueous_component(component)), + solid_concentration(vars.m_solid_var.aqueous_component(component)), + partial_pressure(vars.get_pressure_main_variables(component)), + saturation(vars.m_liquid_var.water()), + porosity(vars.m_porosity), + liquid_diffusivity(vars.m_liquid_diffusivity), + relative_liquid_diffusivity(vars.m_relative_liquid_diffusivity), + advection_flux(vars.m_advection_flux) +{} + +PressureVariableBox::PressureVariableBox(UnsaturatedVariables& vars, index_t component): + partial_pressure(vars.get_pressure_main_variables(component)), + liquid_saturation(vars.m_liquid_var.water()), + porosity(vars.m_porosity), + gas_diffusivity(vars.m_gas_diffusivity), + relative_gas_diffusivity(vars.m_relative_gas_diffusivity) +{} + +// ====================== +// Unsaturated variables +// ====================== + +UnsaturatedVariables::UnsaturatedVariables( + mesh::Mesh1DPtr the_mesh, + database::RawDatabasePtr the_database, + std::vector has_gas + ): + database::DatabaseHolder(the_database), + m_mesh(the_mesh), + m_has_gas(has_gas), + // main var + m_liquid_var(the_database->nb_component(),the_mesh->nb_nodes()), + m_gas_var(the_database->nb_component(),the_mesh->nb_nodes(), has_gas), + m_solid_var(the_database->nb_component(),the_mesh->nb_nodes()), + // chem sol + m_chem_sol(the_mesh->nb_nodes()), + // secondary vars + m_porosity(the_mesh->nb_nodes()), + m_water_aq_concentration(the_mesh->nb_nodes()), + m_liquid_permeability(the_mesh->nb_nodes()), + m_relative_liquid_permeability(the_mesh->nb_nodes()), + m_liquid_diffusivity(the_mesh->nb_nodes()), + m_relative_liquid_diffusivity(the_mesh->nb_nodes()), + m_gas_diffusivity(the_mesh->nb_nodes()), + m_relative_gas_diffusivity(the_mesh->nb_nodes()), + m_capillary_pressure(the_mesh->nb_nodes()), + m_advection_flux(the_mesh->nb_nodes()) +{ + specmicp_assert(has_gas.size() == (unsigned) the_database->nb_component()); + // + // trick to avoid more branching in equations : + m_gas_var.hard_reset(1, the_mesh->nb_nodes()); +} + + +void UnsaturatedVariables::reset_main_variables() +{ + m_liquid_var.reset(); + m_gas_var.reset(); + m_solid_var.reset(); +} + +SaturationVariableBox +UnsaturatedVariables::get_saturation_variables() { + return SaturationVariableBox(*this); +} + +PressureVariableBox +UnsaturatedVariables::get_vapor_pressure_variables() { + return PressureVariableBox(*this, 0); +} + + +LiquidAqueousComponentVariableBox +UnsaturatedVariables::get_liquid_aqueous_component_variables( + index_t component + ) { + return LiquidAqueousComponentVariableBox(*this, component); +} + +PressureVariableBox +UnsaturatedVariables::get_pressure_variables( + index_t component + ) { + return PressureVariableBox(*this, component); +} + + +MainVariable& UnsaturatedVariables::get_pressure_main_variables(index_t component) +{ + if (component_has_gas(component)) + { + specmicp_assert(m_gas_var(component) != nullptr); + return *m_gas_var(component); + } + else + return *m_gas_var(1); // a zero var..., this is a trick to always get a valid info +} + +index_t UnsaturatedVariables::nb_gas() { + return std::count(m_has_gas.begin(), m_has_gas.end(), true); +} + +void UnsaturatedVariables::set_relative_variables( + index_t node + ) +{ + const scalar_t saturation = m_liquid_var.water().variable(node); + + m_relative_liquid_diffusivity(node) = m_relative_liquid_diffusivity_f(node, saturation); + m_relative_gas_diffusivity(node) = m_relative_gas_diffusivity_f(node, saturation); + m_relative_liquid_permeability(node) = m_relative_liquid_permeability_f(node, saturation); + m_capillary_pressure(node) = m_capillary_pressure_f(node, saturation); + +} + +void UnsaturatedVariables::set_relative_variables() +{ + for (index_t node: m_mesh->range_nodes()) + { + set_relative_variables(node); + } +} + + + +} //end namespace unsatura +} //end namespace systems +} //end namespace reactmicp +} //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/variables.hpp b/src/reactmicp/systems/unsaturated/variables.hpp new file mode 100644 index 0000000..c5e40b7 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/variables.hpp @@ -0,0 +1,257 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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_REACTMICP_UNSATURATED_VARIABLES_HPP +#define SPECMICP_REACTMICP_UNSATURATED_VARIABLES_HPP + +#include "../../../types.hpp" +#include "../../../dfpm/meshes/mesh1d.hpp" +#include "../../../database/database_holder.hpp" + +#include "variables_sub.hpp" + +#include "../../solver/staggers_base/variables_base.hpp" + +#include + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + + +// forward declarations +// declared in variables_box.hpp +struct SaturationVariableBox; +struct PressureVariableBox; +struct LiquidAqueousComponentVariableBox; + +//! \brief A function dependent on the saturation +using user_model_saturation_f = std::function; + +class SPECMICP_DLL_PUBLIC UnsaturatedVariables: + public database::DatabaseHolder, + public solver::VariablesBase +{ +public: + UnsaturatedVariables( + mesh::Mesh1DPtr the_mesh, + database::RawDatabasePtr the_database, + std::vector has_gas); + + //! Reset the main variables in case of failure + void reset_main_variables() override; + + //! \brief Return the number of governing equations + //! + //! The governing equations are the transport equations + index_t nb_governing_equations() { + return get_database()->nb_aqueous_components()+1; + } + + //! \brief Return the number of nodes + index_t nb_nodes() { + return m_mesh->nb_nodes(); + } + + //! \brief Return the mesh + mesh::Mesh1DPtr get_mesh() { + return m_mesh; + } + + //! \brief Return the variables for the saturation equation + SaturationVariableBox get_saturation_variables(); + //! \brief Return the variables for the vapor pressure equation + PressureVariableBox get_vapor_pressure_variables(); + //! \brief Return the variables for the liquid transport of an aqueous component + LiquidAqueousComponentVariableBox get_liquid_aqueous_component_variables(index_t component); + //! \brief Return the variables for the gaseous transport of an aqueous component + PressureVariableBox get_pressure_variables(index_t component); + + + //! \brief Set the capillary pressure model + void set_capillary_pressure_model(user_model_saturation_f func) { + m_capillary_pressure_f = func; + } + + //! \brief Set the vapor pressure model + void set_vapor_pressure_model(user_model_saturation_f func) { + m_vapor_pressure_f = func; + } + + //! \brief Set the relative liquid permeability model + void set_relative_liquid_permeability_model(user_model_saturation_f func) { + m_relative_liquid_permeability_f = func; + } + + //! \brief Set the relative liquid diffusivity model + void set_relative_liquid_diffusivity_model(user_model_saturation_f func) { + m_relative_liquid_diffusivity_f = func; + } + + //! \brief Set the relative gas diffusivity model + void set_relative_gas_diffusivity_model(user_model_saturation_f func) { + m_relative_gas_diffusivity_f = func; + } + + //! \brief Set the relative variables at 'node' + void set_relative_variables(index_t node); + + //! \brief Set the relative variables + void set_relative_variables(); + + //! \brief Return a reference to the porosity variables + SecondaryTransientVariable& get_porosity() { + return m_porosity;} + //! \brief Return a reference to the water aq. concentration variables + SecondaryTransientVariable& get_water_aqueous_concentration() { + return m_water_aq_concentration;} + //! \brief Return a reference to the liquid permeability + SecondaryVariable& get_liquid_permeability() { + return m_liquid_permeability; + } + //! \brief Return a reference to the relative liquid permeability + SecondaryVariable& get_relative_permeability() { + return m_relative_liquid_permeability; + } + //! \brief Return a reference to the liquid diffusivity + SecondaryVariable& get_liquid_diffusivity() { + return m_liquid_diffusivity; + } + //! \brief Return a reference to the relative liquid diffusivity + SecondaryVariable& get_relative_liquid_diffusivity() { + return m_relative_liquid_diffusivity; + } + //! \brief Return a reference to the gas diffusivity + SecondaryVariable& get_gas_diffusivity() { + return m_gas_diffusivity; + } + //! \brief Return a reference to the relative gas diffusivity + SecondaryVariable& get_relative_gas_diffusivity() { + return m_relative_gas_diffusivity; + } + SecondaryVariable& get_advection_flux() { + return m_advection_flux; + } + + //! \brief Return the liquid saturation + MainVariable& get_liquid_saturation() { + return m_liquid_var.water(); + } + + //! \brief Return the variables corresponding to the pressure + //! \param component water or aqueous component + MainVariable& get_pressure_main_variables(index_t component); + + //! \brief Return the aqueous concentration of an aqueous component + MainVariable& get_aqueous_concentration(index_t aq_component) { + specmicp_assert(aq_component >= 2 and aq_component < get_database()->nb_component()); + return m_liquid_var.aqueous_component(aq_component); + } + + //! \brief Return the solid variables + MainVariable& get_solid_concentration(index_t component) { + specmicp_assert_component_bounds(component, get_database()); + return m_solid_var.component(component); + } + + //! \brief Return true if a component has a gas equation + //! \param component water or aqueous component + bool component_has_gas(index_t component) { + specmicp_assert_component_bounds(component, get_database()); + return m_has_gas[component]; + } + + //! \brief Return the number of gas in the system + index_t nb_gas(); + +private: + mesh::Mesh1DPtr m_mesh; + + std::vector m_has_gas; + + // main variables + // -------------- + ListMainVariable m_liquid_var; + ListMainVariable m_gas_var; + ListMainVariable m_solid_var; + + // chemistry + // --------- + ChemistrySolutions m_chem_sol; + + // secondary variables + // ------------------- + SecondaryTransientVariable m_porosity; + SecondaryTransientVariable m_water_aq_concentration; + + SecondaryVariable m_liquid_permeability; + SecondaryVariable m_relative_liquid_permeability; + SecondaryVariable m_liquid_diffusivity; + SecondaryVariable m_relative_liquid_diffusivity; + SecondaryVariable m_gas_diffusivity; + SecondaryVariable m_relative_gas_diffusivity; + SecondaryVariable m_capillary_pressure; + + SecondaryVariable m_advection_flux; + + // user models + // ----------- + user_model_saturation_f m_capillary_pressure_f {nullptr}; + user_model_saturation_f m_vapor_pressure_f {nullptr}; + + user_model_saturation_f m_relative_liquid_permeability_f {nullptr}; + user_model_saturation_f m_relative_liquid_diffusivity_f {nullptr}; + user_model_saturation_f m_relative_gas_diffusivity_f {nullptr}; + + // These struct are declared as friends, + // To make their initialisation easier + // (initialisation of referenece, and const reference) + friend struct SaturationVariableBox; + friend struct PressureVariableBox; + friend struct LiquidAqueousComponentVariableBox; +}; + +//! \brief Shared pointer the unsaturated variables +using UnsaturatedVariablesPtr = std::shared_ptr; + +//! \brief Cast from base variables +inline UnsaturatedVariablesPtr cast_from_base(std::shared_ptr ptr) { + return std::static_pointer_cast(ptr); +} + +} //end namespace unsaturated +} //end namespace systems +} //end namespace reactmicp +} //end namespace specmicp + +#endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLES_HPP diff --git a/src/reactmicp/systems/unsaturated/variables_box.hpp b/src/reactmicp/systems/unsaturated/variables_box.hpp new file mode 100644 index 0000000..ef24160 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/variables_box.hpp @@ -0,0 +1,137 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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_REACTMICP_UNSATURATED_VARIABLESBOX +#define SPECMICP_REACTMICP_UNSATURATED_VARIABLESBOX + +#include "../../../types.hpp" +#include "variables_sub.hpp" + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + +//! \brief Base class for a box of variables +//! \internal +struct SPECMICP_DLL_LOCAL BaseVariableBox +{ +}; + +//! \brief Variables box for the saturation equation +//! \internal +struct SPECMICP_DLL_LOCAL SaturationVariableBox: + public BaseVariableBox +{ + MainVariable& liquid_saturation; + const MainVariable& solid_concentration; + const MainVariable& vapor_pressure; + + const SecondaryTransientVariable& aqueous_concentration; + + const SecondaryTransientVariable& porosity; + const SecondaryVariable& liquid_permeability; + const SecondaryVariable& liquid_diffusivity; + + SecondaryVariable& relative_liquid_permeability; + SecondaryVariable& relative_liquid_diffusivity; + SecondaryVariable& capillary_pressure; + + SecondaryVariable& advection_flux; + + user_model_saturation_f capillary_pressure_f; + user_model_saturation_f relative_liquid_permeability_f; + user_model_saturation_f relative_liquid_diffusivity_f; + + //! + //! \brief SaturationVariableBox + //! \param vars The main set of variables + //! + SaturationVariableBox(UnsaturatedVariables& vars); + // implementation fo the constructor is in variables.cpp +}; + +//! \brief Set of variables for the liquid transport equation of an aqueous component +//! \internal +struct SPECMICP_DLL_LOCAL LiquidAqueousComponentVariableBox: + public BaseVariableBox +{ + MainVariable& aqueous_concentration; + const MainVariable& solid_concentration; + const MainVariable& partial_pressure; + + const MainVariable& saturation; + + + const SecondaryTransientVariable& porosity; + const SecondaryVariable& liquid_diffusivity; + const SecondaryVariable& relative_liquid_diffusivity; + const SecondaryVariable& advection_flux; + + //! + //! \brief LiquidAqueousComponentVariableBox + //! \param vars The main set of variables + //! \param component an aqueous component + //! + LiquidAqueousComponentVariableBox(UnsaturatedVariables& vars, index_t component); + // implementation fo the constructor is in variables.cpp +}; + +//! \brief Set of variables for the gas transport equation of an aqueous component +//! \internal +struct SPECMICP_DLL_LOCAL PressureVariableBox: + public BaseVariableBox +{ + MainVariable& partial_pressure; + + const MainVariable& liquid_saturation; + + const SecondaryTransientVariable& porosity; + const SecondaryVariable& gas_diffusivity; + const SecondaryVariable& relative_gas_diffusivity; + + //! + //! \brief PressureVariableBox + //! \param vars The main set of variables + //! \param component The component (can be water) + //! + PressureVariableBox(UnsaturatedVariables& vars, index_t component); + // implementation fo the constructor is in variables.cpp +}; + + +} //end namespace unsaturated +} //end namespace systems +} //end namespace reactmicp +} //end namespace specmicp + +#endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLESBOX diff --git a/src/reactmicp/systems/unsaturated/variables_sub.hpp b/src/reactmicp/systems/unsaturated/variables_sub.hpp new file mode 100644 index 0000000..05de2ee --- /dev/null +++ b/src/reactmicp/systems/unsaturated/variables_sub.hpp @@ -0,0 +1,225 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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 variables_sub.hpp +//! \brief Substruct for the variables + +#ifndef SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP +#define SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP + +#include "../../../types.hpp" +#include "../../../specmicp/adimensional/adimensional_system_solution.hpp" + +#include +#include + +#include "../../../utils/compat.hpp" + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + + +struct SPECMICP_DLL_LOCAL BaseVariable +{ + + Vector variable; //!< The variables + + //! \brief Return value of the variable + scalar_t& operator() (index_t node) { + return variable(node); + } + + //! \brief Return value of the variable + scalar_t operator() (index_t node) const { + return variable(node); + } + + index_t size() const { + return variable.rows(); + } + + void set_constant(scalar_t value) { + variable.setConstant(value); + } + + void set_zero() { + variable.setZero(); + } + + BaseVariable(index_t size): + variable(Vector::Zero(size)) + {} +}; + + +struct SPECMICP_DLL_LOCAL BaseTransientVariable: + public BaseVariable +{ + Vector velocity; //!< Rate of change of the variable + Vector predictor; //!< Predictor, value before first iteration + + BaseTransientVariable(index_t size): + BaseVariable(size), + velocity(Vector::Zero(size)), + predictor(Vector::Zero(size)) + {} +}; + +//! \brief A main variable, used to solve the governing equations +struct SPECMICP_DLL_LOCAL MainVariable: + public BaseTransientVariable +{ + Vector transport_fluxes; //!< Transport fluxes of the corresponding governing equations + Vector chemistry_rate; //!< Chemical exchange rate + + MainVariable(index_t size): + BaseTransientVariable(size), + transport_fluxes(Vector::Zero(size)), + chemistry_rate(Vector::Zero(size)) + {} + + //! \brief Reset the variables + //! + //! This function is called if the solver failed and must be restarted + void reset() { + transport_fluxes.setZero(); + chemistry_rate.setZero(); + velocity.setZero(); + variable = predictor; + } +}; + +//! \brief List of main variables +struct SPECMICP_DLL_LOCAL ListMainVariable +{ + // Variables are stored as pointer so we can skip components + // that do not exist + // But still access their variables through their index + using value_type = std::unique_ptr; + using vector_type = std::vector; + using iterator = vector_type::iterator; + using const_iterator = vector_type::const_iterator; + + vector_type variables; + + MainVariable& water() { + return *variables[0]; + } + + MainVariable& aqueous_component(index_t aq_component) { + specmicp_assert(variables[aq_component] != nullptr); + return *variables[aq_component]; + } + + MainVariable& component(index_t component) { + specmicp_assert(variables[component] != nullptr); + return *variables[component]; + } + + MainVariable* operator() (index_t component) { + return variables[component].get(); + } + + // Implementations of the following methods + // is in variables.cpp + + //! \brief Initialize the variables + ListMainVariable(index_t nb_vars, index_t nb_nodes); + //! \brief Initialize some variables + //! + //! Only initialize values given by 'to_init' + ListMainVariable(index_t nb_vars, index_t nb_nodes, std::vector to_init); + + //! \brief Reset the variables + void reset(); + + //! \brief Hard reset, erase all info for a component; + void hard_reset(index_t component, index_t nb_nodes); +}; + +//! \brief A secondary variable +struct SPECMICP_DLL_LOCAL SecondaryVariable: + public BaseVariable +{ + SecondaryVariable(index_t size): + BaseVariable(size) + {} + +}; + +//! \brief A secondary transient variable +//! +//! e.g. the porosity +struct SPECMICP_DLL_LOCAL SecondaryTransientVariable: + public BaseTransientVariable +{ + SecondaryTransientVariable(index_t size): + BaseTransientVariable(size) + {} +}; + +//! \brief The chemistry solutions +struct SPECMICP_DLL_LOCAL ChemistrySolutions +{ + std::vector solutions; + + ChemistrySolutions(index_t size): + solutions(size) + {} + + //! \brief Return a solution + AdimensionalSystemSolution& solution(index_t node) { + return solutions[node]; + } + + //! \brief Return a const solution + const AdimensionalSystemSolution& solution(index_t node) const { + return solutions[node]; + } + + //! \brief Set a solution + void update_solution(index_t node, const AdimensionalSystemSolution& new_solution) { + solutions[node] = new_solution; + } +}; + + + + +} //end namespace unsaturated +} //end namespace systems +} //end namespace reactmicp +} //end namespace specmicp + +#endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP diff --git a/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt index 8ad6388..1187a3f 100644 --- a/tests/reactmicp/CMakeLists.txt +++ b/tests/reactmicp/CMakeLists.txt @@ -1,44 +1,60 @@ # Mesh # ---- set(test_meshes_files meshes/test_meshes.cpp meshes/test_uniform_mesh_1d.cpp meshes/test_generic1d.cpp meshes/test_axisymmetric_uniform_mesh_1d.cpp meshes/test_configuration.cpp ) add_catch_test(NAME meshes SOURCES "${test_meshes_files}" LINK_LIBRARIES dfpm_static specmicp_common_static ${YAML_LIBRARIES}) FILE(COPY meshes/unif_mesh.yaml meshes/ramp_mesh.yaml meshes/unif_axisym.yaml DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) # Reactmicp : Reactive transport solver # ------------------------------------- set(test_reactmicp_files solver/test_reactive_transport_solver.cpp solver/test_coupling.cpp ) add_catch_test(NAME reactmicp SOURCES ${test_reactmicp_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) # Saturated diffusion using new reactive transport solver # ------------------------------------------------------- set(test_reactmicp_saturated_files systems/saturated_react/test_reactmicp_saturated_react.cpp systems/saturated_react/speciation_system.cpp systems/saturated_react/variables.cpp systems/saturated_react/transport_program.cpp systems/saturated_react/equilibrium_stagger.cpp ) add_catch_test(NAME reactmicp_saturated_react SOURCES ${test_reactmicp_saturated_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) + +# Unsaturated reactive transport system +# ------------------------------------- + +set(test_reactmicp_unsaturated_files + systems/unsaturated/test_reactmicp_unsaturated.cpp + systems/unsaturated/variables.cpp +) + +set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\") +set_source_files_properties(systems/unsaturated/variables.cpp PROPERTIES COMPILE_DEFINITIONS +"TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}" +) + +add_catch_test(NAME reactmicp_unsaturated SOURCES ${test_reactmicp_unsaturated_files} + LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) diff --git a/tests/reactmicp/systems/unsaturated/test_reactmicp_unsaturated.cpp b/tests/reactmicp/systems/unsaturated/test_reactmicp_unsaturated.cpp new file mode 100644 index 0000000..178916e --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/test_reactmicp_unsaturated.cpp @@ -0,0 +1,3 @@ +#define CATCH_CONFIG_MAIN +#include "catch.hpp" + diff --git a/tests/reactmicp/systems/unsaturated/variables.cpp b/tests/reactmicp/systems/unsaturated/variables.cpp new file mode 100644 index 0000000..dd9eb2c --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/variables.cpp @@ -0,0 +1,116 @@ +#include "catch.hpp" + +#include "reactmicp/systems/unsaturated/variables.hpp" +#include "reactmicp/systems/unsaturated/variables_box.hpp" +#include "dfpm/meshes/mesh1d.hpp" +#include "dfpm/meshes/uniform_mesh1d.hpp" +#include "database/database.hpp" + +using namespace specmicp; +using namespace specmicp::reactmicp::systems::unsaturated; + + +static mesh::Mesh1DPtr get_mesh(index_t nb_nodes) +{ + return mesh::uniform_mesh1d(nb_nodes, 1.0, 1.0); +} + + +static specmicp::database::RawDatabasePtr get_database() +{ + static database::RawDatabasePtr raw_data {nullptr}; + if (raw_data == nullptr) + { + specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); + thedatabase.keep_only_components({"H2O", "H[+]", "Ca[2+]", "Si(OH)4"}); + raw_data = thedatabase.get_database(); + raw_data->freeze_db(); + } + return raw_data; +} + +TEST_CASE("Variables Types", "[variables],[types]") { + SECTION("ListMainVariables") { + ListMainVariable alist(5, 4); + + REQUIRE(alist.water().size() == 4); + REQUIRE(alist.aqueous_component(4).size() == 4); + REQUIRE(alist(1) == nullptr); + + + ListMainVariable alist2(4, 6, {true, false, true, false}); + + REQUIRE(alist2.water().size() == 6); + REQUIRE(alist2.aqueous_component(2).size() == 6); + REQUIRE(alist2(1) == nullptr); + REQUIRE(alist2(3) == nullptr); + } +} + +TEST_CASE("Variables", "[variables]") { + + SECTION("Initialisation") { + mesh::Mesh1DPtr the_mesh = get_mesh(10); + database::RawDatabasePtr raw_db = get_database(); + REQUIRE(raw_db->is_valid()); + UnsaturatedVariables vars(the_mesh, raw_db, {true, false, false, false, false}); + + REQUIRE(vars.get_database()->is_valid()); + REQUIRE(vars.get_database()->get_hash() == raw_db->get_hash()); + + REQUIRE(vars.nb_nodes() == 10); + REQUIRE(vars.nb_governing_equations() == 4); + } + + SECTION("Saturation variable") { + index_t nb_nodes = 10; + + mesh::Mesh1DPtr the_mesh = get_mesh(nb_nodes); + database::RawDatabasePtr raw_db = get_database(); + UnsaturatedVariables vars(the_mesh, raw_db, {true, false, false, false, false}); + + SaturationVariableBox satvar = vars.get_saturation_variables(); + + REQUIRE(satvar.liquid_saturation.size() == nb_nodes); + + satvar.liquid_saturation(3) = 2.0; + CHECK(satvar.liquid_saturation(3) == 2.0); + + satvar.liquid_saturation.velocity(2) = 3.0; + CHECK(satvar.liquid_saturation.velocity(2) == 3.0); + + } + + SECTION("Vapor pressure variable") { + index_t nb_nodes = 10; + + mesh::Mesh1DPtr the_mesh = get_mesh(nb_nodes); + database::RawDatabasePtr raw_db = get_database(); + UnsaturatedVariables vars(the_mesh, raw_db, {true, false, false, false, false}); + + PressureVariableBox vap_pressure = vars.get_vapor_pressure_variables(); + REQUIRE(vap_pressure.partial_pressure.size() == nb_nodes); + } + + SECTION("liquid aqueous component variable") { + index_t nb_nodes = 10; + + mesh::Mesh1DPtr the_mesh = get_mesh(nb_nodes); + database::RawDatabasePtr raw_db = get_database(); + UnsaturatedVariables vars(the_mesh, raw_db, {true, false, false, false, false}); + + LiquidAqueousComponentVariableBox liq_var = vars.get_liquid_aqueous_component_variables(2); + REQUIRE(liq_var.aqueous_concentration.size() == nb_nodes); + } + + SECTION("gaseous aqueous component variable") { + index_t nb_nodes = 10; + + mesh::Mesh1DPtr the_mesh = get_mesh(nb_nodes); + database::RawDatabasePtr raw_db = get_database(); + UnsaturatedVariables vars(the_mesh, raw_db, {true, false, false, false, false}); + + PressureVariableBox gas_var = vars.get_pressure_variables(2); + REQUIRE(gas_var.partial_pressure.size() == nb_nodes); // should still work + } +}