PAvgCalculator.hpp

/*
  Copyright 2020, 2023 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 PAVG_CALCULATOR_HPP
#define PAVG_CALCULATOR_HPP

#include <array>
#include <cstddef>
#include <functional>
#include <memory>
#include <optional>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace Opm {

class Connection;
class GridDims;
class PAvg;
template<class Scalar> class PAvgDynamicSourceData;
class WellConnections;

} // namespace Opm

namespace Opm {

template <typename Scalar>
class PAvgCalculatorResult;

template <typename Scalar>
PAvgCalculatorResult<Scalar>
linearCombination(const Scalar alpha, PAvgCalculatorResult<Scalar>        x,
                  const Scalar beta , const PAvgCalculatorResult<Scalar>& y);

template <typename Scalar>
class PAvgCalculatorResult
{
private:
    template <typename T>
    friend PAvgCalculatorResult<T>
    linearCombination(const T alpha, PAvgCalculatorResult<T> x,
                      const T beta , const PAvgCalculatorResult<T>& y);

public:
    enum class WBPMode
    {
        WBP,  //< Connecting cells
        WBP4, //< Immediate neighbours
        WBP5, //< Connecting cells and immediate neighbours
        WBP9, //< Connecting cells, immediate, and diagonal neighbours
    };

    PAvgCalculatorResult& set(const WBPMode type, const Scalar wbp)
    {
        this->wbp_[this->index(type)] = wbp;

        return *this;
    }

    Scalar value(const WBPMode type) const
    {
        return this->wbp_[this->index(type)];
    }

private:
    static constexpr auto NumModes =
        static_cast<std::size_t>(WBPMode::WBP9) + 1;

    using WBPStore = std::array<Scalar, NumModes>;

    WBPStore wbp_{};

    constexpr typename WBPStore::size_type index(const WBPMode mode) const
    {
        return static_cast<typename WBPStore::size_type>(mode);
    }
};

template<class Scalar>
class PAvgCalculator
{
protected:
    class Accumulator;

public:
    class Sources
    {
    public:
        Sources& wellBlocks(const PAvgDynamicSourceData<Scalar>& wbSrc)
        {
            this->wb_ = &wbSrc;
            return *this;
        }

        Sources& wellConns(const PAvgDynamicSourceData<Scalar>& wcSrc)
        {
            this->wc_ = &wcSrc;
            return *this;
        }

        const PAvgDynamicSourceData<Scalar>& wellBlocks() const
        {
            return *this->wb_;
        }

        const PAvgDynamicSourceData<Scalar>& wellConns() const
        {
            return *this->wc_;
        }

    private:
        const PAvgDynamicSourceData<Scalar>* wb_{nullptr};

        const PAvgDynamicSourceData<Scalar>* wc_{nullptr};
    };

    PAvgCalculator(const GridDims&        cellIndexMap,
                   const WellConnections& connections);

    virtual ~PAvgCalculator();

    void pruneInactiveWBPCells(const std::vector<bool>& isActive);

    void inferBlockAveragePressures(const Sources& sources,
                                    const PAvg&    controls,
                                    const Scalar   gravity,
                                    const Scalar   refDepth);

    const std::vector<std::size_t>& allWBPCells() const
    {
        return this->contributingCells_;
    }

    std::vector<std::size_t> allWellConnections() const;

    const PAvgCalculatorResult<Scalar>& averagePressures() const
    {
        return this->averagePressures_;
    }

protected:
    class Accumulator
    {
    public:
        using LocalRunningAverages = std::array<Scalar, 8>;

        Accumulator();

        ~Accumulator();

        Accumulator(const Accumulator& rhs);

        Accumulator(Accumulator&& rhs);

        Accumulator& operator=(const Accumulator& rhs);

        Accumulator& operator=(Accumulator&& rhs);

        Accumulator& addCentre(const Scalar weight,
                               const Scalar press);

        Accumulator& addRectangular(const Scalar weight,
                                    const Scalar press);

        Accumulator& addDiagonal(const Scalar weight,
                                 const Scalar press);

        Accumulator& add(const Scalar       weight,
                         const Accumulator& other);

        void prepareAccumulation();

        void prepareContribution();

        void commitContribution(const Scalar innerWeight = -1.0);

        // Please note that member functions \c getRunningAverages() and \c
        // assignRunningAverages() are concessions to parallel/MPI runs, and
        // especially for simulation runs with distributed wells.  In this
        // situation we need a way to access, communicate/collect/sum, and
        // assign partial results.  Moreover, the \c LocalRunningAverages
        // should be treated opaquely apart from applying a global reduction
        // operation.  In other words, the intended/expected use case is
        //
        //   Accumulator a{}
        //   ...
        //   auto avg = a.getRunningAverages()
        //   MPI_Allreduce(avg, MPI_SUM)
        //   a.assignRunningAverages(avg)
        //
        // Any other use is probably a bug and the above is the canonical
        // implementation of member function collectGlobalContributions() in
        // MPI aware sub classes of PAvgCalculator.

        LocalRunningAverages getRunningAverages() const;

        void assignRunningAverages(const LocalRunningAverages& avg);

        PAvgCalculatorResult<Scalar> getFinalResult() const;

    private:
        class Impl;

        std::unique_ptr<Impl> pImpl_;
    };

    Accumulator accumCTF_{};

    Accumulator accumPV_{};

private:
    using ContrIndexType = std::vector<std::size_t>::size_type;

    using SetupMap = std::unordered_map<std::size_t, ContrIndexType>;

    enum class NeighbourKind
    {
        Rectangular,

        Diagonal,
    };

    struct PAvgConnection
    {
        PAvgConnection(const Scalar         ctf_arg,
                       const ContrIndexType cell_arg)
            : ctf  { ctf_arg }
            , cell { cell_arg }
        {}

        Scalar ctf{};

        ContrIndexType cell{};

        std::vector<ContrIndexType> rectNeighbours{};

        std::vector<ContrIndexType> diagNeighbours{};
    };

    typename std::vector<PAvgConnection>::size_type numInputConns_{};

    std::vector<PAvgConnection> connections_{};

    std::vector<typename std::vector<PAvgConnection>::size_type> openConns_{};

    std::vector<typename std::vector<PAvgConnection>::size_type> inputConn_{};

    std::vector<std::size_t> contributingCells_{};

    PAvgCalculatorResult<Scalar> averagePressures_{};

    void addConnection(const GridDims&   cellIndexMap,
                       const Connection& conn,
                       SetupMap&         setupHelperMap);

    void pruneInactiveConnections(const std::vector<bool>& isActive);

    void accumulateLocalContributions(const Sources& sources,
                                      const PAvg&    controls,
                                      const Scalar   gravity,
                                      const Scalar   refDepth);

    virtual void collectGlobalContributions();

    void assignResults(const PAvg& controls);

    void addNeighbour(std::optional<std::size_t> neighbour,
                      NeighbourKind              neighbourKind,
                      SetupMap&                  setupHelperMap);

    std::size_t lastConnsCell() const;

    void addNeighbours_X(const GridDims& grid, SetupMap& setupHelperMap);

    void addNeighbours_Y(const GridDims& grid, SetupMap& setupHelperMap);

    void addNeighbours_Z(const GridDims& grid, SetupMap& setupHelperMap);

    template <typename ConnIndexMap, typename CTFPressureWeightFunction>
    void accumulateLocalContributions(const Sources&             sources,
                                      const PAvg&                controls,
                                      const Scalar               gravity,
                                      const Scalar               refDepth,
                                      const std::vector<Scalar>& connDensity,
                                      ConnIndexMap               connIndex,
                                      CTFPressureWeightFunction  ctfPressWeight);

    template <typename ConnIndexMap>
    void accumulateLocalContributions(const Sources&             sources,
                                      const PAvg&                controls,
                                      const Scalar               gravity,
                                      const Scalar               refDepth,
                                      const std::vector<Scalar>& connDensity,
                                      ConnIndexMap&&             connIndex);

    void accumulateLocalContribOpen(const Sources&             sources,
                                    const PAvg&                controls,
                                    const Scalar               gravity,
                                    const Scalar               refDepth,
                                    const std::vector<Scalar>& connDP);

    void accumulateLocalContribAll(const Sources&             sources,
                                   const PAvg&                controls,
                                   const Scalar               gravity,
                                   const Scalar               refDepth,
                                   const std::vector<Scalar>& connDensity);

    template <typename ConnIndexMap>
    std::vector<Scalar>
    connectionDensityWell(const std::size_t nconn,
                          const Sources&    sources,
                          ConnIndexMap      connIndex) const;

    template <typename ConnIndexMap>
    std::vector<Scalar>
    connectionDensityRes(const std::size_t nconn,
                         const Sources&    sources,
                         ConnIndexMap      connIndex) const;

    std::vector<Scalar>
    connectionDensity(const Sources& sources,
                      const PAvg&    controls,
                      const Scalar   gravity) const;
};

} // namespace Opm

#endif // PAVG_CALCULATOR_HPP