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 && HAVE_DUNE_ISTL
 
#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