InterRegFlow.hpp

/*
  Copyright (c) 2022 Equinor ASA

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_OUTPUT_DATA_INTERREGFLOW_HPP
#define OPM_OUTPUT_DATA_INTERREGFLOW_HPP

#include <algorithm>
#include <cassert>
#include <array>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <type_traits>
#include <utility>

namespace Opm { namespace data {

    template <typename RandIt>
    class InterRegFlow
    {
    public:
        using ElmT = std::remove_cv_t<
            std::remove_reference_t<
                typename std::iterator_traits<RandIt>::value_type
            >>;

        enum class Component : char {
            Oil, Gas, Water, Disgas, Vapoil,

            // Must be last enumerator
            NumComponents,
        };

        enum class Direction : char {
            Positive, Negative,
        };

        class FlowRates
        {
        public:
            FlowRates()
            {
                this->rate_.fill(ElmT{});
            }

            ElmT& operator[](const Component i)
            {
                return this->rate_[this->index(i)];
            }

            friend class InterRegFlow;

        private:
            std::array<ElmT, static_cast<std::size_t>(Component::NumComponents)> rate_{};

            std::size_t index(const Component i) const
            {
                return static_cast<std::size_t>(i);
            }
        };

        explicit InterRegFlow(RandIt begin, RandIt end)
            : elements_(begin, end)
        {}

        InterRegFlow(const InterRegFlow&) = delete;

        InterRegFlow(InterRegFlow&& rhs)
            : elements_(rhs.elements_.first, rhs.elements_.second)
        {
            rhs.elements_.second = rhs.elements_.first; // rhs -> empty
        }

        InterRegFlow& operator=(const InterRegFlow& rhs)
        {
            this->copyIn(rhs);

            return *this;
        }

        InterRegFlow& operator=(InterRegFlow&& rhs)
        {
            if (! this->isValid()) {
                this->elements_ = rhs.elements_;
            }
            else {
                this->copyIn(rhs);
            }

            rhs.elements_.second = rhs.elements_.first; // rhs -> empty

            return *this;
        }

        template <typename OtherRandIt>
        std::enable_if_t<
            std::is_convertible_v<typename InterRegFlow<OtherRandIt>::ElmT, ElmT>,
        InterRegFlow&> operator+=(const InterRegFlow<OtherRandIt>& rhs)
        {
            std::ranges::transform(*this, rhs, this->begin(), std::plus<>{});

            return *this;
        }

        template <typename OtherRandIt>
        std::enable_if_t<
            !std::is_same_v<RandIt, OtherRandIt> &&
            std::is_convertible_v<typename InterRegFlow<OtherRandIt>::ElmT, ElmT>,
        InterRegFlow&> operator=(const InterRegFlow<OtherRandIt>& rhs)
        {
            this->copyIn(rhs.begin(), rhs.end());

            return *this;
        }

        void addFlow(const ElmT sign, const FlowRates& q)
        {
            assert (this->isValid());

            const auto numComp = static_cast<std::size_t>(Component::NumComponents);

            for (auto component = 0*numComp; component < numComp; ++component) {
                this->add(sign * q.rate_[component], component);
            }
        }

        constexpr static std::size_t bufferSize() noexcept
        {
            return InterRegFlow::index(Component::NumComponents, Direction::Positive);
        }

        constexpr ElmT flow(const Component component) const noexcept
        {
            // Add components since Positive and Negative are stored as
            // signed quantities.  In other words flow(x, Negative) <= 0
            // while flow(x, Positive) >= 0).
            return this->flow(component, Direction::Positive)
                +  this->flow(component, Direction::Negative);
        }

        constexpr ElmT flow(const Component component,
                            const Direction direction) const noexcept
        {
            return *(this->elements_.first + InterRegFlow::index(component, direction));
        }

        constexpr bool empty() const noexcept
        {
            return this->begin() == this->end();
        }

        constexpr bool isValid() const noexcept
        {
            using sz_t = decltype(InterRegFlow::bufferSize());

            const auto& [begin, end] = this->elements_;

            return static_cast<sz_t>(std::distance(begin, end))
                == InterRegFlow::bufferSize();
        }

        RandIt begin() const noexcept
        {
            return this->elements_.first;
        }

        RandIt end() const noexcept
        {
            return this->elements_.second;
        }

    private:
        std::pair<RandIt, RandIt> elements_;

        constexpr static std::size_t
        index(const std::size_t component, const Direction direction)
        {
            return 2*component + (direction == Direction::Negative);
        }

        constexpr static std::size_t
        index(const Component component, const Direction direction)
        {
            return InterRegFlow::index(static_cast<std::size_t>(component), direction);
        }

        void add(const ElmT rate, const std::size_t component)
        {
            const auto direction = std::signbit(rate)
                ? Direction::Negative : Direction::Positive;

            auto* rateVec = &*this->elements_.first;
            rateVec[InterRegFlow::index(component, direction)] += rate;
        }

        void copyIn(const InterRegFlow& rhs)
        {
            if (this->elements_ != rhs.elements_) {
                this->copyIn(rhs.elements_.first, rhs.elements_.second);
            }
        }

        template <typename OtherRandIt>
        void copyIn(OtherRandIt begin, OtherRandIt end)
        {
            std::copy(begin, end, this->elements_.first);
        }
    };

}} // namespace Opm::data

#endif // OPM_OUTPUT_DATA_INTERREGFLOW_HPP