diff --git a/src/reactmicp/systems/unsaturated/variables.cpp b/src/reactmicp/systems/unsaturated/variables.cpp index fce45d3..687ca93 100644 --- a/src/reactmicp/systems/unsaturated/variables.cpp +++ b/src/reactmicp/systems/unsaturated/variables.cpp @@ -1,269 +1,332 @@ -/*------------------------------------------------------------------------------- +/*------------------------------------------------------------------------------ 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" +#include "../../../utils/log.hpp" + +#include + 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), constants(vars.m_constants), 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): +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), constants(vars.m_constants) {} -PressureVariableBox::PressureVariableBox(UnsaturatedVariables& vars, index_t component): +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), constants(vars.m_constants) {} void ConstantBox::scale(units::LengthUnit length_unit) { scalar_t scaling_factor = units::scaling_factor(length_unit); viscosity_liquid_water *= 1.0/scaling_factor; rt *= 1.0/(scaling_factor*scaling_factor); + total_pressure *= 1.0/scaling_factor; } // ====================== // Unsaturated variables // ====================== +namespace internal { + +//! \brief Find the correspond gas id for a component with gas +static index_t get_id_gas(database::RawDatabasePtr the_database, + index_t component) +{ + index_t id = no_species; + for (auto gas: the_database->range_gas()) + { + const scalar_t nu = the_database->nu_gas(gas, component); + if (nu == 0.0) + continue; + + if (nu != 1.0) + { + ERROR << "Database not in the good format" + << "; expect : gas <=> component \n" + << "component : " + << the_database->get_label_component(component) << " \n" + << "gas : " + << the_database->get_label_gas(gas) << "."; + throw std::runtime_error("Badly formatted database"); + } + else + { + id = gas; + } + } + if (id == no_species) + { + ERROR << "No corresponding gas in the database for component : " + << the_database->get_label_component(component); + throw std::runtime_error("Badly formatted database"); + } + return id; +} + +} //end namespace internal + 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), + m_id_gas(the_database->nb_component(), no_species), // 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()); + + for (auto component: the_database->range_aqueous_component()) + { + if (has_gas[component]) + { + m_id_gas[component] = internal::get_id_gas(the_database, component); + } + } } UnsaturatedVariables::UnsaturatedVariables( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, std::vector has_gas, units::LengthUnit length_unit): UnsaturatedVariables(the_mesh, the_database, has_gas) { m_constants.scale(length_unit); } 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); + 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 index 40e0949..c1146cc 100644 --- a/src/reactmicp/systems/unsaturated/variables.hpp +++ b/src/reactmicp/systems/unsaturated/variables.hpp @@ -1,286 +1,366 @@ /*------------------------------------------------------------------------------- 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 +//! \file variables.hpp +//! \name Variables for the unsaturated system + #include "types_fwd.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 +//! \brief The unsaturated system +//! +//! The system solve the transport reactive in unsaturated porous medium namespace unsaturated { // forward declarations // declared in variables_box.hpp struct SaturationVariableBox; struct PressureVariableBox; struct LiquidAqueousComponentVariableBox; +//! \brief The unsaturated variables +//! +//! This class contains all the variables for the unsaturated system +//! +//! Note : it should be initialized using the ::VariablesInterface class. class SPECMICP_DLL_PUBLIC UnsaturatedVariables: public database::DatabaseHolder, public solver::VariablesBase { public: + //! \brief Constructor + //! \internal UnsaturatedVariables( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, std::vector has_gas); + //! \brief Constructor + //! \internal UnsaturatedVariables( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, std::vector has_gas, units::LengthUnit length_unit); - //! Reset the main variables in case of failure + + //! \brief 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 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(); + + // Variables Box + // ============= + + //! \name VariableBox + //! \brief This set of variables are passed to the equations + //! + //! They contain all the variables a governing equation should need + //! @{ + //! \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 + + //! \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 pressure diffusion equation + //! + //! Valid for water vapor or the gas corresponding to an aqueous component PressureVariableBox get_pressure_variables(index_t component); + //! @} + + // User model + // ========== + + //! \name User models + //! \brief Models given by the user to compute saturated-dependant variables + //! + //! A model is a function taking a node and the saturation as the argument + //! @{ + //! \brief Set the capillary pressure model void set_capillary_pressure_model(user_model_saturation_f func) { m_capillary_pressure_f = func; } //! \brief Return the capillary pressure model user_model_saturation_f get_capillary_pressure_model() { return m_capillary_pressure_f; } //! \brief Set the vapor pressure model void set_vapor_pressure_model(user_model_saturation_f func) { m_vapor_pressure_f = func; } //! \brief Return the vapor pressure model user_model_saturation_f get_vapor_pressure_model() { return m_vapor_pressure_f; } //! \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(); + //! @} + + + + // Variables Getter + // ================ + + //! \name Variables Getter + //! \brief Return the 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 capillary pressure SecondaryVariable& get_capillary_pressure() { return m_capillary_pressure; } //! \brief Return a reference to the relative liquid permeability - SecondaryVariable& get_relative_permeability() { + SecondaryVariable& get_relative_liquid_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(); - //! \brief Return the R*T product //! //! R : ideal gas constant //! T : temperature scalar_t get_rt() {return m_constants.rt;} + //! \brief return the total pressure + scalar_t get_total_pressure() {return m_constants.total_pressure;} + + //! \brief return the id of a gas + scalar_t get_id_gas(index_t component) { + specmicp_assert_component_bounds(component, get_database()); + return m_id_gas[component]; + } + + AdimensionalSystemSolution& get_adim_solution(index_t node) { + return m_chem_sol.solution(node); + } + + void set_adim_solution(index_t node, + const AdimensionalSystemSolution& new_solution + ) { + m_chem_sol.update_solution(node, new_solution); + } + + //! @} + private: mesh::Mesh1DPtr m_mesh; std::vector m_has_gas; + std::vector m_id_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; // scaled constants // ---------------- ConstantBox m_constants; // 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) { +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 index 85aee22..9102c6a 100644 --- a/src/reactmicp/systems/unsaturated/variables_box.hpp +++ b/src/reactmicp/systems/unsaturated/variables_box.hpp @@ -1,141 +1,144 @@ /*------------------------------------------------------------------------------- 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 +//! \file unsaturated/variables_box.hpp +//! \brief Set of variables passed to the equations + #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; const ConstantBox& constants; 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; const ConstantBox& constants; //! //! \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; const ConstantBox& constants; //! //! \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_interface.cpp b/src/reactmicp/systems/unsaturated/variables_interface.cpp index e5243de..c6e8d81 100644 --- a/src/reactmicp/systems/unsaturated/variables_interface.cpp +++ b/src/reactmicp/systems/unsaturated/variables_interface.cpp @@ -1,521 +1,540 @@ /*------------------------------------------------------------------------------- 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_interface.hpp" #include "variables.hpp" #include "variables_box.hpp" #include "../../solver/staggers_base/variables_base.hpp" #include "../../../utils/compat.hpp" #include // This is a boring class, with only small differences between // all the functions and variables. But, they make all the // difference at the end. The main purpose of the class is to // hide all the variables structure and only present a unified // interface to the user. // None of the structure defined in "variables_sub.hpp" should // be leaked outside this class namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { // typedef // ------- using VariablesBase = solver::VariablesBase; using VariablesBasePtr = std::shared_ptr; // ================================= // // Declaration of the implementation // // ================================= // This class just store and retrieve variables group struct VariablesInterface::VariablesInterfaceImpl { UnsaturatedVariablesPtr m_vars; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_database; VariablesInterfaceImpl( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, const std::vector& component_with_gas ); MainVariable& liquid_saturation() { return m_vars->get_liquid_saturation(); } const MainVariable& liquid_saturation() const { return m_vars->get_liquid_saturation(); } MainVariable& aqueous_concentration(index_t component) { return m_vars->get_aqueous_concentration(component); } const MainVariable& aqueous_concentration(index_t component) const { return m_vars->get_aqueous_concentration(component); } MainVariable& partial_pressure(index_t component) { if (not m_vars->component_has_gas(component)) { throw std::logic_error("Component has no associated gas"); } return m_vars->get_pressure_main_variables(component); } const MainVariable& partial_pressure(index_t component) const { return m_vars->get_pressure_main_variables(component); } MainVariable& solid_concentration(index_t component) { return m_vars->get_solid_concentration(component); } const MainVariable& solid_concentration(index_t component) const { return m_vars->get_solid_concentration(component); } SecondaryTransientVariable& porosity() { return m_vars->get_porosity(); } const SecondaryTransientVariable& porosity() const { return m_vars->get_porosity(); } SecondaryTransientVariable& water_aqueous_concentration() { return m_vars->get_water_aqueous_concentration(); } const SecondaryTransientVariable& water_aqueous_concentration() const { return m_vars->get_water_aqueous_concentration(); } SecondaryVariable& liquid_diffusivity() { return m_vars->get_liquid_diffusivity(); } const SecondaryVariable& liquid_diffusivity() const { return m_vars->get_liquid_diffusivity(); } SecondaryVariable& liquid_permeability() { return m_vars->get_liquid_permeability(); } const SecondaryVariable& liquid_permeability() const { return m_vars->get_liquid_permeability(); } SecondaryVariable& gas_diffusivity() { return m_vars->get_gas_diffusivity(); } const SecondaryVariable& gas_diffusivity() const { return m_vars->get_gas_diffusivity(); } SecondaryVariable& advection_flux() { return m_vars->get_advection_flux(); } const SecondaryVariable& advection_flux() const { return m_vars->get_advection_flux(); } }; // ================================= // // Implementation // // ================================= VariablesInterface::VariablesInterface( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, const std::vector& component_with_gas): m_impl(make_unique(the_mesh, the_database, component_with_gas)) { } VariablesInterface::VariablesInterfaceImpl::VariablesInterfaceImpl( mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, const std::vector& component_with_gas ): m_mesh(the_mesh), m_database(the_database) { std::vector comp_has_gas(the_database->nb_component(), false); for (auto component: component_with_gas) { comp_has_gas[component] = true; } m_vars = std::make_shared(the_mesh, the_database, comp_has_gas); } VariablesInterface::~VariablesInterface() = default; std::shared_ptr VariablesInterface::get_base_variables() { return std::static_pointer_cast(m_impl->m_vars); } std::shared_ptr VariablesInterface::get_variables() { return m_impl->m_vars; } UnsaturatedVariables* VariablesInterface::get_raw_variables() { return m_impl->m_vars.get(); } // Saturation const Vector& VariablesInterface::get_liquid_saturation() const { return m_impl->liquid_saturation().variable; } void VariablesInterface::set_liquid_saturation(index_t node, scalar_t value) { m_impl->liquid_saturation().variable(node) = value; m_impl->liquid_saturation().predictor(node) = value; } void VariablesInterface::set_liquid_saturation(const Vector& value) { m_impl->liquid_saturation().variable = value; m_impl->liquid_saturation().predictor = value; } void VariablesInterface::set_liquid_saturation(scalar_t value) { m_impl->liquid_saturation().variable.setConstant(value); m_impl->liquid_saturation().predictor.setConstant(value); } // Aqueous component concentration const Vector& VariablesInterface::get_aqueous_concentration(index_t component) const { return m_impl->aqueous_concentration(component).variable; } void VariablesInterface::set_aqueous_concentration(index_t component, index_t node, scalar_t value) { MainVariable& aq_conc = m_impl->aqueous_concentration(component); aq_conc.variable(node) = value; aq_conc.predictor(node) = value; } void VariablesInterface::set_aqueous_concentration(index_t component, const Vector& value) { MainVariable& aq_conc = m_impl->aqueous_concentration(component); aq_conc.variable = value; aq_conc.predictor = value; } void VariablesInterface::set_aqueous_concentration(index_t component, scalar_t value) { MainVariable& aq_conc = m_impl->aqueous_concentration(component); aq_conc.variable.setConstant(value); aq_conc.predictor.setConstant(value); } // Partial pressure const Vector& VariablesInterface::get_partial_pressure(index_t component) const { return m_impl->partial_pressure(component).variable; } void VariablesInterface::set_partial_pressure(index_t component, index_t node, scalar_t value) { MainVariable& comp_pressure = m_impl->partial_pressure(component); comp_pressure.variable(node) = value; comp_pressure.predictor(node) = value; } void VariablesInterface::set_partial_pressure(index_t component, const Vector& value) { MainVariable& comp_pressure = m_impl->partial_pressure(component); comp_pressure.variable = value; comp_pressure.predictor = value; } void VariablesInterface::set_partial_pressure(index_t component, scalar_t value) { MainVariable& comp_pressure = m_impl->partial_pressure(component); comp_pressure.variable.setConstant(value); comp_pressure.predictor.setConstant(value); } // Solid concentration const Vector& VariablesInterface::get_solid_concentration(index_t component) const { return m_impl->solid_concentration(component).variable; } void VariablesInterface::set_solid_concentration(index_t component, index_t node, scalar_t value) { MainVariable& solid_conc = m_impl->solid_concentration(component); solid_conc.variable(node) = value; solid_conc.predictor(node) = value; } void VariablesInterface::set_solid_concentration(index_t component, const Vector& value) { MainVariable& solid_conc = m_impl->solid_concentration(component); solid_conc.variable = value; solid_conc.predictor = value; } void VariablesInterface::set_solid_concentration(index_t component, scalar_t value) { MainVariable& solid_conc = m_impl->solid_concentration(component); solid_conc.variable.setConstant(value); solid_conc.predictor.setConstant(value); } // Porosity const Vector& VariablesInterface::get_porosity() const { return m_impl->porosity().variable; } void VariablesInterface::set_porosity(index_t node, scalar_t value) { m_impl->porosity().variable(node) = value; m_impl->porosity().predictor(node) = value; } void VariablesInterface::set_porosity(const Vector& value) { m_impl->porosity().variable = value; m_impl->porosity().predictor = value; } void VariablesInterface::set_porosity(scalar_t value) { m_impl->porosity().variable.setConstant(value); m_impl->porosity().predictor.setConstant(value); } // Water aqueous concentration const Vector& VariablesInterface::get_water_aqueous_concentration() const { return m_impl->water_aqueous_concentration().variable; } void VariablesInterface::set_water_aqueous_concentration(index_t node, scalar_t value) { m_impl->water_aqueous_concentration().variable(node) = value; m_impl->water_aqueous_concentration().predictor(node) = value; } void VariablesInterface::set_water_aqueous_concentration(const Vector& value) { m_impl->water_aqueous_concentration().variable = value; m_impl->water_aqueous_concentration().predictor = value; } void VariablesInterface::set_water_aqueous_concentration(scalar_t value) { m_impl->water_aqueous_concentration().variable.setConstant(value); m_impl->water_aqueous_concentration().predictor.setConstant(value); } // Liquid diffusion coefficient const Vector& VariablesInterface::get_liquid_diffusivity() const { return m_impl->liquid_diffusivity().variable; } void VariablesInterface::set_liquid_diffusivity(index_t node, scalar_t value) { m_impl->liquid_diffusivity().variable(node) = value; } void VariablesInterface::set_liquid_diffusivity(const Vector& value) { m_impl->liquid_diffusivity().variable = value; } void VariablesInterface::set_liquid_diffusivity(scalar_t value) { m_impl->liquid_diffusivity().variable.setConstant(value); } // Liquid permeability const Vector& VariablesInterface::get_liquid_permeability() const { return m_impl->liquid_permeability().variable; } void VariablesInterface::set_liquid_permeability(index_t node, scalar_t value) { m_impl->liquid_permeability().variable(node) = value; } void VariablesInterface::set_liquid_permeability(const Vector& value) { m_impl->liquid_permeability().variable = value; } void VariablesInterface::set_liquid_permeability(scalar_t value) { m_impl->liquid_permeability().variable.setConstant(value); } // Gas diffusivity const Vector& VariablesInterface::get_gas_diffusivity() const { return m_impl->gas_diffusivity().variable; } void VariablesInterface::set_gas_diffusivity(index_t node, scalar_t value) { m_impl->gas_diffusivity().variable(node) = value; } void VariablesInterface::set_gas_diffusivity(const Vector& value) { m_impl->gas_diffusivity().variable = value; } void VariablesInterface::set_gas_diffusivity(scalar_t value) { m_impl->gas_diffusivity().variable.setConstant(value); } const Vector& VariablesInterface::get_advection_flux() const { return m_impl->advection_flux().variable; } void VariablesInterface::set_advection_flux(index_t node, scalar_t value) { m_impl->advection_flux().variable(node) = value; } void VariablesInterface::set_advection_flux(const Vector& value) { m_impl->advection_flux().variable = value; } void VariablesInterface::set_advection_flux(scalar_t value) { m_impl->advection_flux().variable.setConstant(value); } // User Models +const Vector& VariablesInterface::get_capillary_pressure() const +{ + return m_impl->m_vars->get_capillary_pressure().variable; +} void VariablesInterface::set_capillary_pressure_model( user_model_saturation_f capillary_pressure_model ) { m_impl->m_vars->set_capillary_pressure_model(capillary_pressure_model); const Vector& sat = get_liquid_saturation(); Vector& cap_pressure = m_impl->m_vars->get_capillary_pressure().variable; for (auto node: m_impl->m_mesh->range_nodes()) { cap_pressure(node) = capillary_pressure_model(node, sat(node)); } } +const Vector& VariablesInterface::get_relative_liquid_diffusivity() const +{ + return m_impl->m_vars->get_relative_liquid_diffusivity().variable; +} + void VariablesInterface::set_relative_liquid_diffusivity_model( user_model_saturation_f relative_liquid_diffusivity_model ) { m_impl->m_vars->set_relative_liquid_diffusivity_model(relative_liquid_diffusivity_model); const Vector& sat = get_liquid_saturation(); Vector& rel_diff = m_impl->m_vars->get_relative_liquid_diffusivity().variable; for (auto node: m_impl->m_mesh->range_nodes()) { rel_diff(node) = relative_liquid_diffusivity_model(node, sat(node)); } } +const Vector& VariablesInterface::get_relative_liquid_permeability() const +{ + return m_impl->m_vars->get_relative_liquid_permeability().variable; +} + void VariablesInterface::set_relative_liquid_permeability_model( user_model_saturation_f relative_liquid_permeability_model ) { m_impl->m_vars->set_relative_liquid_permeability_model(relative_liquid_permeability_model); const Vector& sat = get_liquid_saturation(); - Vector& rel_diff = m_impl->m_vars->get_relative_permeability().variable; + Vector& rel_diff = m_impl->m_vars->get_relative_liquid_permeability().variable; for (auto node: m_impl->m_mesh->range_nodes()) { rel_diff(node) = relative_liquid_permeability_model(node, sat(node)); } } +const Vector& VariablesInterface::get_relative_gas_diffusivity() const +{ + return m_impl->m_vars->get_relative_gas_diffusivity().variable; +} + void VariablesInterface::set_relative_gas_diffusivity_model( user_model_saturation_f relative_gas_diffusivity_model ) { m_impl->m_vars->set_relative_gas_diffusivity_model(relative_gas_diffusivity_model); const Vector& sat = get_liquid_saturation(); Vector& rel_diff = m_impl->m_vars->get_relative_gas_diffusivity().variable; for (auto node: m_impl->m_mesh->range_nodes()) { rel_diff(node) = relative_gas_diffusivity_model(node, sat(node)); } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/variables_interface.hpp b/src/reactmicp/systems/unsaturated/variables_interface.hpp index b40f6a6..dce3b89 100644 --- a/src/reactmicp/systems/unsaturated/variables_interface.hpp +++ b/src/reactmicp/systems/unsaturated/variables_interface.hpp @@ -1,196 +1,351 @@ -/*------------------------------------------------------------------------------- +/*------------------------------------------------------------------------------ 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_VARIABLESINTERFACE_HPP #define SPECMICP_REACTMICP_UNSATURATED_VARIABLESINTERFACE_HPP #include "types_fwd.hpp" #include "../../../dfpm/meshes/mesh1dfwd.hpp" #include "../../../database/database_fwd.hpp" #include //! \file variables_interface.hpp //! \brief Interface to initialize and store variables namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { class UnsaturatedVariables; //! \brief Initialize safely a set of variables class SPECMICP_DLL_PUBLIC VariablesInterface { public: VariablesInterface(mesh::Mesh1DPtr the_mesh, database::RawDatabasePtr the_database, const std::vector& component_with_gas); ~VariablesInterface(); //! \brief Return the variables to pass to the reactive transport solver solver::VariablesBasePtr get_base_variables(); //! \brief Return the variables UnsaturatedVariablesPtr get_variables(); //! \brief Return the variables UnsaturatedVariables* get_raw_variables(); // Variables // Main Variables // ============== // Saturation // ---------- + //! \name Liquid saturation + //! \brief The saturation of the liquid phase + //! @{ + + //! \brief Return the liquid saturation const Vector& get_liquid_saturation() const; + //! \brief Set the liquid saturation at 'node' void set_liquid_saturation(index_t node, scalar_t value); - void set_liquid_saturation(const Vector& value); + //! \brief Copy 'values' to the liquid saturation vector + //! + //! \param values a column vector containing the values, + //! it's size is the number of nodes + void set_liquid_saturation(const Vector& values); + //! \brief Set the saturation to be constant accross the mesh void set_liquid_saturation(scalar_t value); + //! @} // Aqueous concentration // --------------------- + //! \name Aqueous concentration + //! \brief Concentration (mol/volume) of the aqueous component + //! @{ + + //! \brief Return the aqueous concentration of 'component' const Vector& get_aqueous_concentration(index_t component) const; - void set_aqueous_concentration(index_t component, index_t node, scalar_t value); - void set_aqueous_concentration(index_t component, const Vector& value); + //! \brief Set the aqueous concentration of 'component' at 'node' + void set_aqueous_concentration( + index_t component, + index_t node, + scalar_t value + ); + //! \brief Copy 'values' into the aqueous concentration vector + //! of 'component' + //! + //! \param values a column vector containing the values, + //! it's size is the number of nodes + void set_aqueous_concentration(index_t component, const Vector& values); + //! \brief Set a constant aqueous concentration for 'component' void set_aqueous_concentration(index_t component, scalar_t value); - // Vapor pressure + //! @} + + // Partial pressure // -------------- + //! \name Partial pressure + //! \brief Partial pressure for a gas + //! @{ + + //! \brief Return the partial pressure for 'component' const Vector& get_partial_pressure(index_t component) const; + //! \brief Set the partial pressure of 'component' at 'node' void set_partial_pressure(index_t component, index_t node, scalar_t value); - void set_partial_pressure(index_t component, const Vector& value); + //! \brief Copy 'values' into the 'component' partial pressure vector + void set_partial_pressure(index_t component, const Vector& values); + //! \brief Set a constant partial pressure for 'component' void set_partial_pressure(index_t component, scalar_t value); + // @} + // Solid concentration // ------------------- + //! \name Solid concentration + //! \brief Solid concentration (moles/volume) + //! @{ + + //! \brief Return the solid concentration of 'component' const Vector& get_solid_concentration(index_t component) const; + //! \brief Set the solid concentration of 'component' at 'node' void set_solid_concentration(index_t component, index_t node, scalar_t value); - void set_solid_concentration(index_t component, const Vector& value); + //! \brief Copy 'values' to the 'component' solid concentration vector + void set_solid_concentration(index_t component, const Vector& values); + //! \brief Set a constant solid concentration for 'component' void set_solid_concentration(index_t component, scalar_t value); + //! @} // Chemistry // --------- // Secondary Variables // =================== // Porosity // -------- + //! \name Porosity + //! \brief The porosity of the porous media + //! @{ + + //! \brief Return the porosity vector const Vector& get_porosity() const; + //! \brief Set the porosity at 'node' void set_porosity(index_t node, scalar_t value); - void set_porosity(const Vector& value); + //! \brief Copy 'values' into the porosity vector + void set_porosity(const Vector& values); + //! \brief Set a constant porosity accross the mesh void set_porosity(scalar_t value); - // Aqueous water concentration + //! @} + + // Water aqueous concentration // --------------------------- + //! \name Water aqueous concentration + //! \brief Aqueous concentration (moles/volume) for the water + + //! \brief Return the water aqueous concentration const Vector& get_water_aqueous_concentration() const; + //! \brief Set the water aqueous concentration at 'node' void set_water_aqueous_concentration(index_t node, scalar_t value); - void set_water_aqueous_concentration(const Vector& value); + //! \brief Copy values into the water aqueous concentration vector + void set_water_aqueous_concentration(const Vector& values); + //! \brief Set the water aqueous concentration to be constant accross + //! the mesh void set_water_aqueous_concentration(scalar_t value); // Liquid diffusivity // ------------------ + //! \name Liquid diffusion coefficient + //! \brief Intrinsic diffusion coefficient in the liquid phase + //! + //! This value only depends on the microstructure, not the saturation + //! @{ + + //! \brief Return the liquid diffusion coefficient const Vector& get_liquid_diffusivity() const; + //! \brief Set the liquid diffusion coefficient at 'node' void set_liquid_diffusivity(index_t node, scalar_t value); - void set_liquid_diffusivity(const Vector& value); + //! \brief Copy 'values' into the liquid diffusion coefficient vector + void set_liquid_diffusivity(const Vector& values); + //! \brief Set the liquid diffusion coefficient to be constant accross + //! the mesh void set_liquid_diffusivity(scalar_t value); + //! @} + // Liquid permeability // ------------------- + //! \name Liquid permeability + //! \brief Intrinsic permeability of the liquid phase + //! + //! This value only depends on the microstructure, not the saturation + //! @{ + + //! \brief Return the liquid permeability const Vector& get_liquid_permeability() const; + //! \brief Set the liquid permeability at 'node' void set_liquid_permeability(index_t node, scalar_t value); - void set_liquid_permeability(const Vector& value); + //! \brief Copy 'values' into the liquid permeability vector + void set_liquid_permeability(const Vector& values); + //! \brief Set the liquid permeability to be constant accross the mesh void set_liquid_permeability(scalar_t value); + //! @} // Gas diffusivity // ------------------ + //! \name Gas diffusion coefficient + //! \brief Intrinsic diffusion coefficient of the gas phase + //! + //! This value only depends on the microstructure, not the saturation + //! @{ + + //! \brief Return the gas diffusion coefficient vector const Vector& get_gas_diffusivity() const; + //! \brief Set the gas diffusion coefficient at 'node' void set_gas_diffusivity(index_t node, scalar_t value); - void set_gas_diffusivity(const Vector& value); + //! \brief Copy 'values' into the gas diffusion coefficient vector + void set_gas_diffusivity(const Vector& values); + //! \brief Set the gas diffusion coefficient to be constant accross the mesh void set_gas_diffusivity(scalar_t value); + //! @} + // Advection flux // --------------- + //! \name Liquid phase velocity + //! \brief The liquid phase velocity due to the capillary pressure + //! @{ + + //! \brief Return the liquid phase velocity vector const Vector& get_advection_flux() const; + //! \brief Set the liquid phase velocity at 'node' void set_advection_flux(index_t node, scalar_t value); + //! \brief Copy 'values' into the liquid phase velocity vector void set_advection_flux(const Vector& value); + //! \brief Set the liquid phase velocity to be constant accross the mesh void set_advection_flux(scalar_t value); + //! @} + // User models // ----------- + //! \name Capillary pressure + //! \brief The macroscopic capillary pressure + //! @{ + + //! \brief Return the capillary pressure values + const Vector& get_capillary_pressure() const; + //! \brief Set the capillary pressure model void set_capillary_pressure_model( user_model_saturation_f capillary_pressure_model ); + //! @} + + //! \name Relative liquid diffusivity + //! \brief The relative diffusion coefficient of the liquid phase + //! @{ + + //! \brief Return the relative liquid diffusion coefficient values + const Vector& get_relative_liquid_diffusivity() const; + //! \brief Set the relative liquid diffusion coefficient model void set_relative_liquid_diffusivity_model( user_model_saturation_f relative_liquid_diffusivity_model ); + //! @} + + //! \name Relative liquid permeability + //! \brief The relative permeability of the liquid phase + //! @{ + + //! \brief Return the relative liquid permeability values + const Vector& get_relative_liquid_permeability() const; + //! \brief Set the relative liquid permeability model void set_relative_liquid_permeability_model( user_model_saturation_f relative_liquid_permeability_model ); + //! @} + + //! \name Relative gas diffusivity + //! \brief The relative diffusion coefficient of the gas phase + //! + //! @{ + + //! \brief Return the relative gas diffusion coefficient values + const Vector& get_relative_gas_diffusivity() const; + //! \brief Set the relative gas diffusion coefficient model void set_relative_gas_diffusivity_model( user_model_saturation_f relative_gas_diffusivity_model ); + + //! @} + private: - struct VariablesInterfaceImpl; // the implementation + //! \brief The implementation class + //! \internal + struct VariablesInterfaceImpl; + //! \brief The implementation details + //! \internal std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLESINTERFACE_HPP diff --git a/src/reactmicp/systems/unsaturated/variables_sub.hpp b/src/reactmicp/systems/unsaturated/variables_sub.hpp index 8d6e4e4..cfc28e0 100644 --- a/src/reactmicp/systems/unsaturated/variables_sub.hpp +++ b/src/reactmicp/systems/unsaturated/variables_sub.hpp @@ -1,243 +1,267 @@ /*------------------------------------------------------------------------------- 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 +//! \file unsaturated/variables_sub.hpp +//! \brief base structures 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 "../../../utils/compat.hpp" #include "../../../physics/constants.hpp" #include "../../../physics/units.hpp" #include #include namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { - -struct SPECMICP_DLL_LOCAL BaseVariable +//! \brief A simple variable +struct BaseVariable { - Vector variable; //!< The variables + Vector variable; //!< The variable //! \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); } + //! \brief Return the size of the vector index_t size() const { return variable.rows(); } + //! \brief The variables are all set to 'value' void set_constant(scalar_t value) { variable.setConstant(value); } + //! \brief Set the variables to zero void set_zero() { variable.setZero(); } + //! \brief Build a variable of size 'size' + //! + //! The varaibles are initialized to zero BaseVariable(index_t size): variable(Vector::Zero(size)) {} }; - -struct SPECMICP_DLL_LOCAL BaseTransientVariable: +//! \brief A transient variable +//! +//! These type of variable store the rate of change of the variable (velocity) +//! and the value at the beginning of the timestep (predictor) +struct 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 Update the variable + void update(scalar_t dt) { + variable = predictor + dt*velocity; + } }; //! \brief A main variable, used to solve the governing equations -struct SPECMICP_DLL_LOCAL MainVariable: +//! +//! These variables also store the chemistry and transport rates for the +//! coupling +struct MainVariable: public BaseTransientVariable { - Vector transport_fluxes; //!< Transport fluxes of the corresponding governing equations - Vector chemistry_rate; //!< Chemical exchange rate + //! \brief Transport fluxes of the corresponding governing equations + Vector transport_fluxes; + //! \brief Chemical exchange rate + Vector chemistry_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 +struct 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: +struct SecondaryVariable: public BaseVariable { SecondaryVariable(index_t size): BaseVariable(size) {} }; //! \brief A secondary transient variable //! //! e.g. the porosity -struct SPECMICP_DLL_LOCAL SecondaryTransientVariable: +struct SecondaryTransientVariable: public BaseTransientVariable { SecondaryTransientVariable(index_t size): BaseTransientVariable(size) {} }; //! \brief The chemistry solutions -struct SPECMICP_DLL_LOCAL ChemistrySolutions +struct 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) { + void update_solution( + index_t node, + const AdimensionalSystemSolution& new_solution + ) { solutions[node] = new_solution; } }; //! \brief List of constants used in the problems -struct SPECMICP_DLL_LOCAL ConstantBox +struct ConstantBox { //! \brief R*T //! //! R : ideal gas constant //! T : temperature scalar_t rt {constants::gas_constant*units::celsius(25.0)}; //! \brief Viscosity of liquid water scalar_t viscosity_liquid_water {constants::water_viscosity}; + //! \brief Total pressure + scalar_t total_pressure {1.01325e5}; //! \brief Scale the constants void scale(units::LengthUnit length_unit); }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_VARIABLESSUB_HPP