diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index c68f822..0cd1341 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,166 +1,167 @@ # ReactMiCP # ============= set(REACTMICP_SOLVER_DIR solver) set(REACTMICP_SYSTEMS_DIR systems) set(REACTMICP_EQUILIBRIUMCURVE equilibrium_curve) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/types_fwd.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_sub.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_box.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/transport_constraints.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ${REACTMICP_SOLVER_DIR}/timestepper.cpp ${REACTMICP_SOLVER_DIR}/runner.cpp # ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/saturation_equation.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/pressure_equation.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/aqueous_equation.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_interface.cpp + ${REACTMICP_SYSTEMS_DIR}/unsaturated/equilibrium_constraints.cpp io/reactive_transport.cpp io/saturated_react.cpp io/configuration.cpp ) set_visibility_flag(${REACTMICP_LIBFILES}) set_pgo_flag(${REACTMICP_LIBFILES}) if(OPENMP_FOUND) if(SPECMICP_USE_OPENMP) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) endif() endif() add_library(objreactmicp OBJECT ${REACTMICP_LIBFILES}) set_property(TARGET objreactmicp PROPERTY POSITION_INDEPENDENT_CODE 1) add_library(reactmicp SHARED $) target_link_libraries(reactmicp dfpm) install(TARGETS reactmicp LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # include files # ------------- set(REACTMICP_STAGGERBASE_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ) set(REACTMICP_SOLVER_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.hpp ${REACTMICP_SOLVER_DIR}/timestepper.hpp ${REACTMICP_SOLVER_DIR}/runner.hpp ) set(REACTMICP_SATURATEDREACT_INCLUDE_LIST ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.hpp ) set(REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp ) set(REACTMICP_IO_INCLUDE_LIST io/reactive_transport.hpp io/saturated_react.hpp io/configuration.hpp ) install(FILES ${REACTMICP_STAGGERBASE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver/staggers_base ) install(FILES ${REACTMICP_SOLVER_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver ) install(FILES ${REACTMICP_SATURATEDREACT_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/saturated_react ) install(FILES ${REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/equilibrium_curve ) install(FILES ${REACTMICP_IO_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/io ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(reactmicp_static STATIC $) install(TARGETS reactmicp_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(reactmicp_static EXCLUDE_FROM_ALL STATIC $) endif() set_target_properties(reactmicp_static PROPERTIES OUTPUT_NAME reactmicp) target_link_libraries(reactmicp_static dfpm_static) diff --git a/src/reactmicp/systems/unsaturated/equilibrium_constraints.cpp b/src/reactmicp/systems/unsaturated/equilibrium_constraints.cpp new file mode 100644 index 0000000..f6ee94d --- /dev/null +++ b/src/reactmicp/systems/unsaturated/equilibrium_constraints.cpp @@ -0,0 +1,216 @@ +/*------------------------------------------------------------------------------ + +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 "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) +{ + index_t nb_nodes = m_impl->m_node_constraints.size(); + specmicp_assert(constraints_for_each_node.size() + == nb_nodes); + for (index_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 < 0 or id_constraint > m_impl->last()) { + ERROR << "Constraints #" << id_constraint << " does not exist."; + throw std::runtime_error("Invalid constraint request"); + } + for (index_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 new file mode 100644 index 0000000..27f74e2 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/equilibrium_constraints.hpp @@ -0,0 +1,142 @@ +/*------------------------------------------------------------------------------ + +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_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