ParallelIstlInformation.hpp

/*
  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
  Copyright 2014, 2015 Statoil ASA
  Copyright 2015 NTNU

  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_PARALLELISTLINFORMATION_HEADER_INCLUDED
#define OPM_PARALLELISTLINFORMATION_HEADER_INCLUDED

#include <vector>

#if HAVE_MPI

#include <algorithm>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>

#include <dune/istl/owneroverlapcopy.hh>

#include <opm/simulators/utils/ParallelCommunication.hpp>

namespace Opm
{

class ParallelISTLInformation
{
public:
    using ParallelIndexSet = Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet;
    using RemoteIndices = Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices;

    ParallelISTLInformation();

    explicit ParallelISTLInformation(MPI_Comm communicator);

    ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
                            const std::shared_ptr<RemoteIndices>& remoteIndices,
                            MPI_Comm communicator);

    ParallelISTLInformation(const ParallelISTLInformation& other);

    std::shared_ptr<ParallelIndexSet> indexSet() const
    {
        return indexSet_;
    }

    std::shared_ptr<RemoteIndices> remoteIndices() const
    {
        return remoteIndices_;
    }

    Parallel::Communication communicator() const
    {
        return communicator_;
    }
    void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices,
                      std::size_t local_component_size = 0,
                      std::size_t num_components = 1) const;

    template<class T>
    void copyOwnerToAll (const T& source, T& dest) const;

    template<class T>
    const std::vector<double>& updateOwnerMask(const T& container) const;

    const std::vector<double>& getOwnerMask() const
    {
        return ownerMask_;
    }

    template<typename Container, typename BinaryOperator, typename T>
    void computeReduction(const Container& container, BinaryOperator binaryOperator,
                          T& value) const;

private:
    template<typename... Containers, typename... BinaryOperators, typename... ReturnValues>
    void computeTupleReduction(const std::tuple<Containers...>& containers,
                               std::tuple<BinaryOperators...>& operators,
                               std::tuple<ReturnValues...>& values) const;

    std::shared_ptr<ParallelIndexSet> indexSet_;
    std::shared_ptr<RemoteIndices> remoteIndices_;
    Parallel::Communication communicator_;
    mutable std::vector<double> ownerMask_;
};

    namespace Reduction
    {
    // the reduction operation.
    template<typename BinaryOperator>
    struct MaskIDOperator
    {
        template<class T, class T1>
        T operator()(const T& t1, const T& t2, const T1& mask)
        {
            return b_(t1, maskValue(t2, mask));
        }
        template<class T, class T1>
        T maskValue(const T& t, const T1& mask)
        {
            return t*mask;
        }
        BinaryOperator& localOperator()
        {
            return b_;
        }
        template<class T>
        T getInitialValue()
        {
            return T{};
        }
    private:
        BinaryOperator b_;
    };

    template<class T>
    struct InnerProductFunctor
    {
        template<class T1>
        T operator()(const T& t1, const T& t2, const T1& mask)
        {
            T masked =  maskValue(t2, mask);
            return t1 + masked * masked;
        }
        template<class T1>
        T maskValue(const T& t, const T1& mask)
        {
            return t*mask;
        }
        std::plus<T> localOperator()
        {
            return std::plus<T>();
        }
        template<class V, std::enable_if_t<std::is_convertible_v<T, V>, int> = 0>
        V getInitialValue()
        {
            return T{};
        }
    };

    // the reduction operation.
    template<typename BinaryOperator>
    struct MaskToMinOperator
    {
        explicit MaskToMinOperator(BinaryOperator b)
        : b_(b)
        {}
        template<class T, class T1>
        T operator()(const T& t1, const T& t2, const T1& mask)
        {
            return b_(t1, maskValue(t2, mask));
        }
        template<class T, class T1>
        T maskValue(const T& t, const T1& mask)
        {
            if( mask )
            {
                return t;
            }
            else
            {
                return getInitialValue<T>();
            }
        }
        template<class T>
        T getInitialValue()
        {
            return std::numeric_limits<T>::lowest();
        }
        BinaryOperator& localOperator()
        {
            return b_;
        }
    private:
        BinaryOperator b_;
    };

    template<typename BinaryOperator>
    struct MaskToMaxOperator
    {
        explicit MaskToMaxOperator(BinaryOperator b)
        : b_(b)
        {}
        template<class T, class T1>
        T operator()(const T& t1, const T& t2, const T1& mask)
        {
            return b_(t1, maskValue(t2, mask));
        }
        template<class T, class T1>
        T maskValue(const T& t, const T1& mask)
        {
            if( mask )
            {
                return t;
            }
            else
            {
                return std::numeric_limits<T>::max();
            }
        }
        BinaryOperator& localOperator()
        {
            return b_;
        }
        template<class T>
        T getInitialValue()
        {
            return std::numeric_limits<T>::max();
        }
    private:
        BinaryOperator b_;
    };
    template<class T>
    MaskIDOperator<std::plus<T> >
    makeGlobalSumFunctor()
    {
        return MaskIDOperator<std::plus<T> >();
    }
    template<class T>
    auto makeGlobalMaxFunctor()
    {
        struct MaxOp
        {
            using result_type = T;
            const result_type& operator()(const T& t1, const T& t2)
            {
                return std::max(t1, t2);
            }
        };
        return MaskToMinOperator(MaxOp());
    }

    namespace detail
    {
        template<typename T, typename Enable = void>
        struct MaxAbsFunctor
        {
            using result_type = T;
            result_type operator()(const T& t1,
                                   const T& t2)
            {
                return std::max(std::abs(t1), std::abs(t2));
            }
        };

        // Specialization for unsigned integers. They need their own
        // version since abs(x) is ambiguous (as well as somewhat
        // meaningless).
        template<typename T>
        struct MaxAbsFunctor<T, typename std::enable_if<std::is_unsigned<T>::value>::type>
        {
            using result_type = T;
            result_type operator()(const T& t1,
                                   const T& t2)
            {
                return std::max(t1, t2);
            }
        };
    }

    template<class T>
    MaskIDOperator<detail::MaxAbsFunctor<T> >
    makeLInfinityNormFunctor()
    {
        return MaskIDOperator<detail::MaxAbsFunctor<T> >();
    }
    template<class T>
    auto
    makeGlobalMinFunctor()
    {
        struct MinOp
        {
            using result_type = T;
            const result_type& operator()(const T& t1, const T& t2)
            {
                return std::min(t1, t2);
            }
        };
        return MaskToMaxOperator(MinOp());
    }
    template<class T>
    InnerProductFunctor<T>
    makeInnerProductFunctor()
    {
        return InnerProductFunctor<T>();
    }
    } // end namespace Reduction
} // end namespace Opm

#endif

namespace Opm
{
template<class T1>
auto
accumulateMaskedValues(const T1& container, const std::vector<double>* maskContainer)
    -> decltype(container[0]*(*maskContainer)[0]);

} // end namespace Opm

#endif