diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index 2acb19a..5071aeb 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,213 +1,211 @@ # ReactMiCP # ============= set(REACTMICP_SOLVER_DIR solver) set(REACTMICP_SYSTEMS_DIR systems) set(REACTMICP_EQUILIBRIUMCURVE equilibrium_curve) set(REACTMICP_IO io) set(REACTMICP_UNSAT_DIR ${REACTMICP_SYSTEMS_DIR}/unsaturated) 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_UNSAT_DIR}/types_fwd.hpp ${REACTMICP_UNSAT_DIR}/variables_sub.hpp ${REACTMICP_UNSAT_DIR}/variables_box.hpp ${REACTMICP_UNSAT_DIR}/transport_constraints.hpp ${REACTMICP_UNSAT_DIR}/fv_1dof_equation.hpp ${REACTMICP_UNSAT_DIR}/fv_1dof_equation.inl ${REACTMICP_UNSAT_DIR}/init_variables.hpp ${REACTMICP_UNSAT_DIR}/upscaling_stagger_factory.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_UNSAT_DIR}/variables.cpp ${REACTMICP_UNSAT_DIR}/saturation_equation.cpp ${REACTMICP_UNSAT_DIR}/saturation_pressure_equation.cpp ${REACTMICP_UNSAT_DIR}/pressure_equation.cpp ${REACTMICP_UNSAT_DIR}/aqueous_equation.cpp ${REACTMICP_UNSAT_DIR}/aqueous_pressure_equation.cpp ${REACTMICP_UNSAT_DIR}/transport_stagger.cpp ${REACTMICP_UNSAT_DIR}/equilibrium_stagger.cpp ${REACTMICP_UNSAT_DIR}/variables_interface.cpp - ${REACTMICP_UNSAT_DIR}/equilibrium_constraints.cpp ${REACTMICP_UNSAT_DIR}/boundary_conditions.cpp ${REACTMICP_IO}/reactive_transport.cpp ${REACTMICP_IO}/saturated_react.cpp ${REACTMICP_IO}/configuration.cpp ) # hdf5 # ----- if (HDF5_FOUND) list(APPEND REACTMICP_LIBFILES ${REACTMICP_IO}/hdf5_unsaturated.cpp) include_directories(${HDF5_INCLUDE_DIRS}) set_source_files_properties( ${REACTMICP_IO}/hdf5_unsaturated.cpp PROPERTIES COMPILE_DEFINITIONS HDF5_DEFINITIONS ) endif() 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" ) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/unsaturated/transport_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/unsaturated/equilibrium_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_UNSATURATED_INCLUDE_LIST ${REACTMICP_UNSAT_DIR}/types_fwd.hpp ${REACTMICP_UNSAT_DIR}/variables_sub.hpp ${REACTMICP_UNSAT_DIR}/variables_box.hpp ${REACTMICP_UNSAT_DIR}/transport_constraints.hpp ${REACTMICP_UNSAT_DIR}/init_variables.hpp ${REACTMICP_UNSAT_DIR}/variables.hpp ${REACTMICP_UNSAT_DIR}/transport_stagger.hpp ${REACTMICP_UNSAT_DIR}/equilibrium_stagger.hpp ${REACTMICP_UNSAT_DIR}/variables_interface.hpp - ${REACTMICP_UNSAT_DIR}/equilibrium_constraints.hpp ${REACTMICP_UNSAT_DIR}/upscaling_stagger_factory.hpp ${REACTMICP_UNSAT_DIR}/boundary_conditions.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_UNSATURATED_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/unsaturated ) 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/equilibrium_constraints.cpp b/src/reactmicp/systems/unsaturated/equilibrium_constraints.cpp deleted file mode 100644 index 1aff632..0000000 --- a/src/reactmicp/systems/unsaturated/equilibrium_constraints.cpp +++ /dev/null @@ -1,218 +0,0 @@ -/* ============================================================================= - - Copyright (c) 2014 - 2016 - F. Georget Princeton University - All rights reserved. - - -Redistribution and use in source and binary forms, with or without modification, -are permitted provided that the following conditions are met: - - 1. Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - 3. Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived from this - software without specific prior written permission. - -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR -ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON -ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * - -============================================================================= */ - -#include "equilibrium_constraints.hpp" - -#include "../../../specmicp/adimensional/adimensional_system_solver_structs.hpp" - -#include "../../../utils/compat.hpp" -#include "../../../utils/log.hpp" - -#include "../../../dfpm/meshes/mesh1d.hpp" - -namespace specmicp { -namespace reactmicp { -namespace systems { -namespace unsaturated { - -struct EquilibriumConstraints::EquilibriumConstraintsImpl -{ - - EquilibriumConstraintsImpl(index_t nb_nodes): - m_node_constraints(nb_nodes, no_equation), - m_is_fixed_node(nb_nodes, false) - { - m_constraints.reserve(3); - } - - EquilibriumConstraintsImpl(index_t nb_nodes, - const AdimensionalSystemSolverOptions& opts): - m_node_constraints(nb_nodes, no_equation), - m_is_fixed_node(nb_nodes, false), - m_options(opts) - { - m_constraints.reserve(3); - } - - size_t last() {return m_constraints.size() - 1;} - - AdimensionalSystemSolverOptions& get_options() {return m_options;} - - std::vector m_constraints; - std::vector m_node_constraints; - - std::vector m_is_fixed_node; - - void set_fixed_node(index_t node) { - m_is_fixed_node[node] = true; - } - bool is_fixed_node(index_t node) { - return m_is_fixed_node[node]; - } - - index_t nb_nodes() {return m_is_fixed_node.size();} - - AdimensionalSystemSolverOptions m_options; -}; - -EquilibriumConstraints::EquilibriumConstraints(index_t nb_nodes): - m_impl(make_unique(nb_nodes)) -{} - - -EquilibriumConstraints::EquilibriumConstraints( - index_t nb_nodes, - const AdimensionalSystemSolverOptions& opts): - m_impl(make_unique(nb_nodes, opts)) -{} - -EquilibriumConstraints::~EquilibriumConstraints() = default; - -std::size_t EquilibriumConstraints::add_constraints( - const AdimensionalSystemConstraints& constraints) -{ - m_impl->m_constraints.push_back(constraints); - return m_impl->last(); -} - - -std::size_t EquilibriumConstraints::new_constraints() -{ - m_impl->m_constraints.emplace_back(); - return m_impl->last(); -} - -AdimensionalSystemConstraints& -EquilibriumConstraints::operator[] (std::size_t index_constraint) -{ - return m_impl->m_constraints[index_constraint]; -} - - -void EquilibriumConstraints::set_constraints( - std::size_t id_constraint, - const std::vector& nodes - ) -{ - specmicp_assert(id_constraint <= m_impl->last()); - for (auto node: nodes) { - m_impl->m_node_constraints[node] = id_constraint; - } -} - -void EquilibriumConstraints::set_constraints( - const std::vector& constraints_for_each_node) -{ - uindex_t nb_nodes = m_impl->m_node_constraints.size(); - specmicp_assert(constraints_for_each_node.size() - == nb_nodes); - for (uindex_t node=0; node m_impl->last()) { - ERROR << "Node " << node << " : " - << "constraints #" << id << " does not exist."; - throw std::runtime_error("Invalid constraint request"); - } - m_impl->m_node_constraints[node] = constraints_for_each_node[node]; - } -} - - -void EquilibriumConstraints::set_constraint_for_all(std::size_t id_constraint) -{ - if (id_constraint > m_impl->last()) { - ERROR << "Constraints #" << id_constraint << " does not exist."; - throw std::runtime_error("Invalid constraint request"); - } - for (uindex_t node=0; node < m_impl->m_node_constraints.size(); ++node) - { - m_impl->m_node_constraints[node] = id_constraint; - } -} - -AdimensionalSystemConstraints& -EquilibriumConstraints::get_constraints(index_t node) -{ - auto id = m_impl->m_node_constraints[node]; - if (id == no_equation) - { - // check that the node has been initialized before doing anything - ERROR << "Node " << node << " : " - << "request for a non initialized constraints !"; - throw std::runtime_error("Invalid constraint request."); - } - return m_impl->m_constraints[id]; -} - -//! \brief Set the options -AdimensionalSystemSolverOptions& EquilibriumConstraints::get_options() -{ - return m_impl->get_options(); -} - -//! \brief Add a fixed node -void EquilibriumConstraints::add_fixed_node(index_t node) -{ - specmicp_assert(node < m_impl->nb_nodes()); - m_impl->set_fixed_node(node); -} - -//! \brief Add a list of fixed nodes -void EquilibriumConstraints::add_fixed_nodes(std::vector fixed_nodes) -{ - for (auto node: fixed_nodes) { - specmicp_assert(node < m_impl->nb_nodes()); - m_impl->set_fixed_node(node); - } -} - -//! \brief Return true if the node is fixed -bool EquilibriumConstraints::is_fixed_node(index_t node) -{ - return m_impl->is_fixed_node(node); -} - -std::shared_ptr equilibrium_constraints( - mesh::Mesh1DPtr the_mesh - ) -{ - index_t nb_nodes = the_mesh->nb_nodes(); - return std::make_shared(nb_nodes); -} - -} //end namespace unsaturated -} //end namespace systems -} //end namespace reactmicp -} //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/equilibrium_constraints.hpp b/src/reactmicp/systems/unsaturated/equilibrium_constraints.hpp deleted file mode 100644 index b324a23..0000000 --- a/src/reactmicp/systems/unsaturated/equilibrium_constraints.hpp +++ /dev/null @@ -1,144 +0,0 @@ -/* ============================================================================= - - Copyright (c) 2014 - 2016 - F. Georget Princeton University - All rights reserved. - - -Redistribution and use in source and binary forms, with or without modification, -are permitted provided that the following conditions are met: - - 1. Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - 3. Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived from this - software without specific prior written permission. - -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR -ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON -ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * - -============================================================================= */ - -#ifndef SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRUMCONSTRAINTS_HPP -#define SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRUMCONSTRAINTS_HPP - -#include "../../../types.hpp" - -#include "../../../dfpm/meshes/mesh1dfwd.hpp" - -#include - -namespace specmicp { - -struct AdimensionalSystemConstraints; -struct AdimensionalSystemSolverOptions; - -namespace reactmicp { -namespace systems { -namespace unsaturated { - -//! \brief Constraints for the equilibrium stagger -class SPECMICP_DLL_PUBLIC EquilibriumConstraints -{ -public: - //! \brief Build a new set of constraints - EquilibriumConstraints(index_t nb_nodes); - //! \brief Build a new set of constraints - EquilibriumConstraints(index_t nb_nodes, - const AdimensionalSystemSolverOptions& opts); - - ~EquilibriumConstraints(); - - //! \brief Add a new constraints to the set of constraints - //! - //! \param constraints A set of constraints - //! \return the index of the set of constraints - //! - //! \sa new_constraints - std::size_t add_constraints( - const AdimensionalSystemConstraints& constraints); - - //! \brief Add a new default constraint to the system - //! - //! \return the index of the set of constraints - //! - //! The constraints can be obtained through operator[] - //! - //! Ex : - //! - //! EquilibriumConstraints constraints(10); - //! auto id = constraints.new_constraints(); - //! auto& constraint_to_modify = constraints[]; - //! // to stuff to constraint_to_modify - //! - std::size_t new_constraints(); - - //! \brief Return the constraint with index 'index_constraint' - //! - //! \param index_constraint is the index return by add_constraints, - //! or new_constraints - AdimensionalSystemConstraints& operator[] (std::size_t index_constraint); - - //! \brief Set the constraint to be 'id_constraint' for each node in 'nodes' - void set_constraints( - std::size_t id_constraint, - const std::vector& nodes - ); - - //! \brief Set constraints for each node - //! - //! \param constraints_for_each_node Vector containing the index of the - //! constraint for each node - void set_constraints( - const std::vector& constraints_for_each_node); - - //! \brief Set the constraint for all node - void set_constraint_for_all(std::size_t id_constraint); - - //! \brief Return the constraint for node (node) - AdimensionalSystemConstraints& get_constraints(index_t node); - - //! \brief Add a fixed node - void add_fixed_node(index_t node); - - //! \brief Add a list of fixed nodes - void add_fixed_nodes(std::vector fixed_nodes); - - //! \brief Return true if the node is fixed - bool is_fixed_node(index_t node); - - //! \brief Set the options - AdimensionalSystemSolverOptions& get_options(); - -private: - struct EquilibriumConstraintsImpl; - std::unique_ptr m_impl; -}; - -//! \brief Return a constraint container for the EquilibriumStagger -std::shared_ptr equilibrium_constraints( - mesh::Mesh1DPtr the_mesh - ); - - - - -} //end namespace unsaturated -} //end namespace systems -} //end namespace reactmicp -} //end namespace specmicp - -#endif // SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRUMCONSTRAINTS_HPP diff --git a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp index 6e10efe..379af8d 100644 --- a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp @@ -1,506 +1,508 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "equilibrium_stagger.hpp" #include "variables.hpp" -#include "equilibrium_constraints.hpp" - #include "../../solver/staggers_base/stagger_structs.hpp" #include "../../../dfpm/meshes/mesh1d.hpp" #include "../../../database/data_container.hpp" #include "../../../utils/compat.hpp" #include "../../../specmicp/adimensional/adimensional_system_solution.hpp" #include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "../../../specmicp/adimensional/adimensional_system_solver.hpp" #include "../../../utils/log.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { // ============================================ // // Declaration of implementation details // // ============================================ struct EquilibriumStagger::EquilibriumStaggerImpl { + using OptionsVector = EquilibriumStagger::OptionsVector; + scalar_t m_dt {-1}; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_database; - std::shared_ptr m_constraints; + std::shared_ptr m_bcs; + std::shared_ptr m_options; EquilibriumStaggerImpl( std::shared_ptr variables, - std::shared_ptr constraints + std::shared_ptr boundary_conditions, + std::shared_ptr options ): m_mesh(variables->get_mesh()), m_database(variables->get_database()), - m_constraints(constraints) + m_bcs(boundary_conditions), + m_options(options) {} scalar_t dt() {return m_dt;} int compute_one_node(index_t node, UnsaturatedVariables * const vars); void compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ); void analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ); void initialize_timestep_one_node(index_t node, UnsaturatedVariables* vars); void initialize(UnsaturatedVariables* vars); void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars); StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); units::UnitsSet& get_units() { - return m_constraints->get_options().units_set; + return m_options->get("default").units_set; } }; using TrueConstPtr = UnsaturatedVariables * const; inline TrueConstPtr cast_to_var(solver::VariablesBase * const var) { return static_cast(var); } // ============================================ // // Implementation // // ============================================ EquilibriumStagger::EquilibriumStagger( std::shared_ptr variables, - std::shared_ptr constraints + std::shared_ptr boundary_conditions, + std::shared_ptr options ): - m_impl(utils::make_pimpl(variables, constraints)) + m_impl(utils::make_pimpl( + variables, boundary_conditions, options)) { } EquilibriumStagger::~EquilibriumStagger() = default; //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize(true_vars); } //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize_timestep(dt, true_vars); } //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables EquilibriumStagger::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); return m_impl->restart_timestep(true_vars); } int EquilibriumStagger::EquilibriumStaggerImpl::compute_one_node( index_t node, UnsaturatedVariables * const vars ) { - AdimensionalSystemConstraints constraints = m_constraints->get_constraints(node); + AdimensionalSystemConstraints constraints = m_bcs->get_constraint(node); compute_total_concentrations(node, constraints.total_concentrations, vars); // set partial pressure model { user_model_saturation_f callable = vars->get_vapor_pressure_model(); water_partial_pressure_f pv_model = [callable, node](scalar_t sat){ return callable(node, sat); }; constraints.set_water_partial_pressure_model(pv_model); } // We use the current value of the saturation to avoid converging // to the previous solution // This is necessary because the update to the saturation is very small const AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node); Vector x; AdimensionalSystemSolver the_solver; if (likely(previous_solution.is_valid)) // should always be true { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, previous_solution, - m_constraints->get_options() + m_options->operator [](node) ); x = previous_solution.main_variables; // we use the new value of the saturation to avoid a change // to small to be detected by the speciaion solver x(0) = porosity*saturation; // x(0) = volume fraction water } else { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, - m_constraints->get_options() + m_options->operator [](node) ); the_solver.initialise_variables(x, porosity*saturation, -3); } // micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { ERROR << "Failed to solve equilibrium at node " << node << ". \n" << "Return code : " << (int) perf.return_code; return -1; } AdimensionalSystemSolution solution = the_solver.get_raw_solution(x); analyse_solution(node, solution, vars); vars->set_adim_solution(node, solution); return 0; } void EquilibriumStagger::EquilibriumStaggerImpl::compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ) { total_concentrations.resize(m_database->nb_component()); const scalar_t saturation = vars->get_liquid_saturation()(node); const scalar_t porosity = vars->get_porosity()(node); const scalar_t rt = vars->get_rt(); // water { const scalar_t ctilde_w = vars->get_water_aqueous_concentration()(node); const scalar_t cbar_w = vars->get_solid_concentration(0)(node); scalar_t c_w = saturation*porosity*ctilde_w + cbar_w; if (vars->component_has_gas(0)) { const scalar_t pv_w = vars->get_pressure_main_variables(0)(node); c_w += (1.0-saturation)*porosity*pv_w/rt; } total_concentrations(0) = c_w; } total_concentrations(1) = 0.0; // aqueous components for (index_t component: m_database->range_aqueous_component()) { const scalar_t ctilde_i = vars->get_aqueous_concentration(component)(node); const scalar_t cbar_i = vars->get_solid_concentration(component)(node); scalar_t c_i = saturation*porosity*ctilde_i + cbar_i; if (vars->component_has_gas(component)) { const scalar_t pv_i = vars->get_pressure_main_variables(component)(node); c_i += (1.0-saturation)*porosity*pv_i/rt; } total_concentrations(component) = c_i; } } void EquilibriumStagger::EquilibriumStaggerImpl::analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ) { AdimensionalSystemSolutionExtractor extr(solution, m_database, get_units()); const scalar_t saturation = extr.saturation_water(); const scalar_t sat_g = 1.0 - saturation; const scalar_t rho_l = extr.density_water(); const scalar_t porosity = extr.porosity(); const scalar_t rt = vars->get_rt(); // porosity SecondaryTransientVariable& porosity_vars = vars->get_porosity(); const scalar_t porosity_0 = porosity_vars.predictor(node); const scalar_t porosity_vel = (porosity - porosity_0)/m_dt; porosity_vars.variable(node) = porosity; porosity_vars.velocity(node) = porosity_vel; // water // saturation MainVariable& satvars = vars->get_liquid_saturation(); const scalar_t saturation_0 = satvars.predictor(node); const scalar_t sat_vel = (saturation - saturation_0)/m_dt; satvars.variable(node) = saturation; satvars.velocity(node) = sat_vel; { // total aqueous concentration SecondaryTransientVariable& cwtilde_vars = vars->get_water_aqueous_concentration(); const scalar_t cwtilde = rho_l*extr.total_aqueous_concentration(0); const scalar_t cwtilde_0 = cwtilde_vars.predictor(node); const scalar_t cwtilde_vel = (cwtilde - cwtilde_0)/m_dt; cwtilde_vars.variable(node) = cwtilde; cwtilde_vars.velocity(node) = cwtilde_vel; // total solid concentration MainVariable& solid_vars = vars->get_solid_concentration(0); const scalar_t cwbar = extr.total_solid_concentration(0); const scalar_t cwbar_0 = solid_vars.predictor(node); const scalar_t cwbar_vel = (cwbar - cwbar_0)/m_dt; solid_vars.variable(node) = cwbar; solid_vars.velocity(node) = cwbar_vel; solid_vars.chemistry_rate(node) = - cwbar_vel; // vapor pressure if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); const scalar_t pv = vars->get_vapor_pressure_model()(node, saturation); const scalar_t pv_0 = pres_vars.predictor(node); const scalar_t pv_vel = (pv - pv_0)/m_dt; pres_vars.variable(node) = pv; pres_vars.velocity(node) = pv_vel; const scalar_t transient = ( (pv*porosity*sat_g) - (pv_0*porosity_0*(1.0-saturation_0)) )/ (rt*m_dt); const scalar_t pv_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pv_chem_rate; } } // aqueous components for (index_t component: m_database->range_aqueous_component()) { // total aqueous concentration MainVariable& aqconc = vars->get_aqueous_concentration(component); const scalar_t ctildei = rho_l*extr.total_aqueous_concentration(component); const scalar_t ctildei_0 = aqconc.predictor(node); const scalar_t ctildei_vel = (ctildei - ctildei_0) / m_dt; aqconc.variable(node) = ctildei; aqconc.velocity(node) = ctildei_vel; // total solid concentration MainVariable& solconc = vars->get_solid_concentration(component); const scalar_t cbari = extr.total_solid_concentration(component); const scalar_t cbari_0 = solconc.predictor(node); const scalar_t cbari_vel = (cbari - cbari_0) / m_dt; solconc.variable(node) = cbari; solconc.velocity(node) = cbari_vel; solconc.chemistry_rate(node) = - cbari_vel; // partial pressure if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); index_t id_g = vars->get_id_gas(component); const scalar_t pi = extr.fugacity_gas(id_g)*vars->get_total_pressure(); const scalar_t pi_0 = pres_vars.predictor(node); const scalar_t pi_vel = (pi - pi_0) / m_dt; pres_vars.variable(node) = pi; pres_vars.velocity(node) = pi_vel; const scalar_t transient = ( (pi*porosity*sat_g) - (pi_0*porosity_0*(1.0-saturation_0)) ) / (rt * m_dt); const scalar_t pi_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pi_chem_rate; } } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep_one_node( index_t node, UnsaturatedVariables * const vars ) { } EquilibriumStagger::StaggerReturnCode EquilibriumStagger::EquilibriumStaggerImpl::restart_timestep( UnsaturatedVariables* vars ) { int sum_retcode = 0; #if SPECMICP_USE_OPENMP #pragma omp parallel for \ reduction(+: sum_retcode) #endif for (index_t node=0; nodenb_nodes(); ++node) { - if (m_constraints->is_fixed_node(node)) continue; + if (not m_bcs->is_normal_node(node)) continue; sum_retcode += compute_one_node(node, vars); } if (sum_retcode != 0) { return StaggerReturnCode::UnknownError; } return StaggerReturnCode::ResidualMinimized; } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables* vars ) { m_dt = dt; // initialize - - - //water { MainVariable& cbar_w = vars->get_solid_concentration(0); cbar_w.predictor = cbar_w.variable; cbar_w.velocity.setZero(); cbar_w.chemistry_rate.setZero(); SecondaryTransientVariable& ctilde_w = vars->get_water_aqueous_concentration(); ctilde_w.predictor = ctilde_w.variable; ctilde_w.velocity.setZero(); if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // aqueous component for (index_t component: m_database->range_aqueous_component()) { MainVariable& cbar_i = vars->get_solid_concentration(component); cbar_i.predictor = cbar_i.variable; cbar_i.velocity.setZero(); cbar_i.chemistry_rate.setZero(); if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // porosity { SecondaryTransientVariable& porosity = vars->get_porosity(); porosity.predictor = porosity.variable; porosity.velocity.setZero(); } for (auto node: m_mesh->range_nodes()) { initialize_timestep_one_node(node, vars); } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize( UnsaturatedVariables * const vars ) { // do nothing by default } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/equilibrium_stagger.hpp b/src/reactmicp/systems/unsaturated/equilibrium_stagger.hpp index 58e57c7..9495ae8 100644 --- a/src/reactmicp/systems/unsaturated/equilibrium_stagger.hpp +++ b/src/reactmicp/systems/unsaturated/equilibrium_stagger.hpp @@ -1,101 +1,118 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRIUMSTAGGER_HPP #define SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRIUMSTAGGER_HPP #include "../../solver/staggers_base/chemistry_stagger_base.hpp" #include "../../../utils/pimpl_ptr.hpp" +#include "boundary_conditions.hpp" +#include "utils/cached_vector.hpp" + namespace specmicp { + +struct AdimensionalSystemSolverOptions; + namespace reactmicp { namespace systems { namespace unsaturated { -class EquilibriumConstraints; class UnsaturatedVariables; class SPECMICP_DLL_PUBLIC EquilibriumStagger: public solver::ChemistryStaggerBase { public: using StaggerReturnCode = solver::StaggerReturnCode; using VariablesBase = solver::VariablesBase; + using OptionsVector = utils::NameCachedVector; EquilibriumStagger( std::shared_ptr variables, - std::shared_ptr constraints + std::shared_ptr boundary_conditions, + std::shared_ptr options ); ~EquilibriumStagger(); + + //! \brief Build a default SolverOptions vector + static std::shared_ptr + make_default_options(index_t nb_nodes, std::string name="default") { + return std::make_shared(nb_nodes, name); + } + //! \brief Return an equilibrium stagger + static std::shared_ptr make( + std::shared_ptr variables, + std::shared_ptr boundary_conditions, + std::shared_ptr options + ) + { + return std::make_shared( + variables,boundary_conditions, options); + } + //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables void initialize(VariablesBase * const var); //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables void initialize_timestep(scalar_t dt, VariablesBase * const var); //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables StaggerReturnCode restart_timestep(VariablesBase * const var); private: struct SPECMICP_DLL_LOCAL EquilibriumStaggerImpl; utils::pimpl_ptr m_impl; }; -//! \brief Return an equilibrium stagger -inline std::shared_ptr equilibrium_stagger( - std::shared_ptr variables, - std::shared_ptr constraints - ) -{ - return std::make_shared(variables, constraints); -} + } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_EQUILIBRIUMSTAGGER_HPP diff --git a/src/reactmicp_unsaturated.hpp b/src/reactmicp_unsaturated.hpp index 479c501..f9750c5 100644 --- a/src/reactmicp_unsaturated.hpp +++ b/src/reactmicp_unsaturated.hpp @@ -1,59 +1,57 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ //! \file reactmicp_unsaturated.hpp //! \brief The unsaturated system of the ReactMiCP solver #ifndef SPECMICP_REACMICP_UNSATURATED_HPP #define SPECMICP_REACMICP_UNSATURATED_HPP // Core files #include "reactmicp_core.hpp" // System #include "reactmicp/systems/unsaturated/variables_interface.hpp" #include "reactmicp/systems/unsaturated/init_variables.hpp" #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" -#include "reactmicp/systems/unsaturated/equilibrium_constraints.hpp" #include "reactmicp/systems/unsaturated/equilibrium_stagger.hpp" - #include "reactmicp/systems/unsaturated/transport_stagger.hpp" // Output #endif // SPECMICP_REACMICP_UNSATURATED_HPP diff --git a/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp b/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp index 08972ff..cdb20d2 100644 --- a/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp +++ b/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp @@ -1,214 +1,193 @@ #include "catch.hpp" #include "database/database.hpp" #include "database/data_container.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "reactmicp/systems/unsaturated/equilibrium_constraints.hpp" #include "reactmicp/systems/unsaturated/equilibrium_stagger.hpp" #include "reactmicp/systems/unsaturated/variables_interface.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" +#include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "utils/log.hpp" +#include "utils/cached_vector.hpp" +#include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; 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", "HCO3[-]"} ); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } static AdimensionalSystemSolution get_init_solution(database::RawDatabasePtr the_database) { Formulation formulation; formulation.set_mass_solution(500); formulation.add_mineral("C2S", 1000.0); formulation.add_mineral("C3S", 4000.0); formulation.add_aqueous_species("CO2", 20); Dissolver dissolver(the_database); Vector tot_conc = dissolver.dissolve(formulation, false); AdimensionalSystemConstraints constraints; constraints.total_concentrations = tot_conc; AdimensionalSystemSolver adim_solver(the_database, constraints); adim_solver.get_options().solver_options.set_tolerance(1e-14); Vector x; adim_solver.initialise_variables(x, 0.8, -3); micpsolver::MiCPPerformance perf = adim_solver.solve(x); specmicp_assert(perf.return_code > micpsolver::MiCPSolverReturnCode::Success); return adim_solver.get_raw_solution(x); } TEST_CASE("EquilibriumConstraints", "[Equilibrium],[Constraints]") { init_logger(&std::cout, logger::Warning); SECTION("Initialisation") { - EquilibriumConstraints constraints(10); - - auto id_0 = constraints.new_constraints(); - CHECK(id_0 == 0); - - AdimensionalSystemConstraints& constraint = constraints[id_0]; - constraint.set_charge_keeper(2); - - AdimensionalSystemConstraints new_constraint; - new_constraint.set_charge_keeper(3); - - auto id_1 = constraints.add_constraints(new_constraint); - CHECK(id_1 == 1); - - CHECK(constraints[id_1].charge_keeper == new_constraint.charge_keeper); - - constraints.set_constraints(id_0, {0, 1, 2, 3, 4}); - constraints.set_constraints(id_1, {5, 6, 7, 8}); - - CHECK(constraints.get_constraints(0).charge_keeper == 2); - CHECK(constraints.get_constraints(1).charge_keeper == 2); - CHECK(constraints.get_constraints(2).charge_keeper == 2); - CHECK(constraints.get_constraints(3).charge_keeper == 2); - CHECK(constraints.get_constraints(4).charge_keeper == 2); - - CHECK(constraints.get_constraints(5).charge_keeper == 3); - CHECK(constraints.get_constraints(6).charge_keeper == 3); - CHECK(constraints.get_constraints(7).charge_keeper == 3); - CHECK(constraints.get_constraints(8).charge_keeper == 3); - - std::cout << "This error is ok : "; - CHECK_THROWS(constraints.get_constraints(9)); - - std::vector id_constraints(10, 1); - id_constraints[0] = 0; - - constraints.set_constraints(id_constraints); - - CHECK(constraints.get_constraints(0).charge_keeper == 2); - CHECK(constraints.get_constraints(1).charge_keeper == 3); - CHECK(constraints.get_constraints(2).charge_keeper == 3); - CHECK(constraints.get_constraints(3).charge_keeper == 3); - CHECK(constraints.get_constraints(4).charge_keeper == 3); - CHECK(constraints.get_constraints(5).charge_keeper == 3); - CHECK(constraints.get_constraints(6).charge_keeper == 3); - CHECK(constraints.get_constraints(7).charge_keeper == 3); - CHECK(constraints.get_constraints(8).charge_keeper == 3); - CHECK(constraints.get_constraints(9).charge_keeper == 3); - - id_constraints[5] = 2; - std::cout << "This error is ok : "; - CHECK_THROWS(constraints.set_constraints(id_constraints)); + auto bcs = BoundaryConditions::make(10); + bcs->add_fixed_node(0); + bcs->add_fixed_node(9); + + AdimensionalSystemConstraints& def_constraint = bcs->get_constraint("default"); + def_constraint.set_charge_keeper(2); + + AdimensionalSystemConstraints& fixed_constraint = + bcs->fork_constraints(0, "fixed_compo"); + fixed_constraint.set_charge_keeper(3); + + + bcs->set_constraint(9, "fixed_compo"); + + CHECK(bcs->get_constraint(0).charge_keeper == 3); + CHECK(bcs->get_constraint(1).charge_keeper == 2); + CHECK(bcs->get_constraint(2).charge_keeper == 2); + CHECK(bcs->get_constraint(3).charge_keeper == 2); + CHECK(bcs->get_constraint(4).charge_keeper == 2); + CHECK(bcs->get_constraint(5).charge_keeper == 2); + CHECK(bcs->get_constraint(6).charge_keeper == 2); + CHECK(bcs->get_constraint(7).charge_keeper == 2); + CHECK(bcs->get_constraint(8).charge_keeper == 2); + CHECK(bcs->get_constraint(9).charge_keeper == 3); } } TEST_CASE("Equilibrium stagger", "[Equilibrium],[Stagger]") { init_logger(&std::cout, logger::Warning); index_t nb_nodes {10}; scalar_t dx {1.0}; scalar_t cross_section {5.0}; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section}); database::RawDatabasePtr raw_data = get_database(); - std::shared_ptr constraints = equilibrium_constraints(the_mesh); - constraints->set_constraint_for_all(constraints->new_constraints()); - constraints->get_options().solver_options.set_tolerance(1e-14); + + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes()); + + + auto options = EquilibriumStagger::make_default_options(the_mesh->nb_nodes()); + options->get("default").solver_options.set_tolerance(1e-14); + index_t id_co2 = raw_data->get_id_component_from_element("C"); VariablesInterface interface(the_mesh, raw_data, {id_co2}); units::UnitsSet units_set; AdimensionalSystemSolution solution = get_init_solution(raw_data); AdimensionalSystemSolutionExtractor extr(solution, raw_data, units_set); scalar_t rho_l = extr.density_water(); interface.set_porosity(extr.porosity()); scalar_t init_saturation = extr.saturation_water(); interface.set_liquid_saturation(init_saturation); scalar_t init_water_aqconc = rho_l*extr.total_aqueous_concentration(0); interface.set_water_aqueous_concentration(init_water_aqconc); scalar_t init_water_solidconc = extr.total_solid_concentration(0); interface.set_solid_concentration(0, init_water_solidconc); for (index_t comp: raw_data->range_aqueous_component()) { interface.set_aqueous_concentration(comp, rho_l*extr.total_aqueous_concentration(comp)); interface.set_solid_concentration(comp, extr.total_solid_concentration(comp)); } // pv co2 UnsaturatedVariablesPtr vars = interface.get_variables(); scalar_t pv_co2 = extr.fugacity_gas(vars->get_id_gas(id_co2) )*vars->get_total_pressure(); std::cout << "pv co2 : " << pv_co2 << " tot pressure : " << vars->get_total_pressure() << std::endl; interface.set_partial_pressure(id_co2, pv_co2); scalar_t init_co2_aqconc = rho_l*extr.total_aqueous_concentration(id_co2); scalar_t init_co2_solidconc = extr.total_solid_concentration(id_co2); for (index_t node: the_mesh->range_nodes()) { vars->set_adim_solution(node, solution); } SECTION("Initialisation") { - EquilibriumStagger stagger(vars, constraints); + EquilibriumStagger stagger(vars, bcs, options); // Check that nothing has changed after we solved an already solved problem stagger.initialize(vars.get()); stagger.initialize_timestep(1.0, vars.get()); solver::StaggerReturnCode retcode = stagger.restart_timestep(vars.get()); CHECK(retcode > solver::StaggerReturnCode::NotConvergedYet); for (index_t node: the_mesh->range_nodes()) { CHECK(interface.get_liquid_saturation()(node) == init_saturation); CHECK(interface.get_water_aqueous_concentration()(node)/init_water_aqconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_solid_concentration(0)(node) / init_water_solidconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_aqueous_concentration(id_co2)(node)/init_co2_aqconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_solid_concentration(id_co2)(node) / init_co2_solidconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_partial_pressure(id_co2)(node) == Approx(pv_co2).epsilon(1e-12)); } } }