average.h

Go to the documentation of this file.

#pragma once
 
#include "prolongation.h"
#ifdef MPI_VERSION
#include "mpi_prolongation.h"
#endif
 
namespace dg{
 
template<class IMatrix, class ContainerType>
struct Average
{
    using container_type = ContainerType;
    using value_type = dg::get_value_type<ContainerType>;
    Average() = default;
 
    template<class Topology, size_t Md>
    Average( const Topology& g, std::array<unsigned, Md> axes)
    {
        m_prolongation = dg::create::prolongation( g, axes);
        m_average      = dg::create::projection( axes, g);
 
        auto remains = dg::create::detail::complement( g, axes);
 
        auto tmp = dg::evaluate( dg::one, g.grid(remains[0]));
        for( unsigned u=1; u<remains.size(); u++)
            tmp = dg::kronecker( dg::Product(), tmp,
                dg::evaluate(dg::one, g.grid(remains[u])));  //note tmp comes first (important for comms in MPI)
        m_tmp = dg::construct<ContainerType>(tmp);
        m_area_inv = 1.;
        for( unsigned u=0; u<Md; u++)
            m_area_inv /= g.l(axes[u]);
    }
 
    template< class Topology, typename = std::enable_if_t<Topology::ndim() == 2>>
    Average( const Topology& g, enum coo2d axes)
    : Average( g, coo2axis(axes))
    {
    }
 
    template< class Topology, typename = std::enable_if_t<Topology::ndim() == 3>>
    Average( const Topology& g, enum coo3d axes)
    : Average( g, coo2axis(axes))
    {
    }
 
    template< class Topology, typename = std::enable_if_t<Topology::ndim() == 3>>
    Average( const Topology& g, std::array<enum coo3d,2> axes)
    : Average( g, {coo2axis(axes[0]), coo2axis(axes[1])})
    {
    }
 
    void operator() (const ContainerType& src, ContainerType& res, bool extend
            = true)
    {
        // Optimization note: Notice tmp is smaller than src so
        // blas1 functions do not matter so much performance wise
        // Remove possible NaN in tmp (because symv does not do it)
        // (This must be done to remove potential NaN from previous computations,
        // otherwise the class can be "infected" with NaN)
        dg::blas1::copy( 0., m_tmp);
        // Cusparse makes no guarantee on NaN behaviour
        // m_average in MPI has the "allreduce" distribution!
        dg::blas2::symv( m_average, src, m_tmp);
        dg::blas1::scal( m_tmp, m_area_inv);
        if( extend )
        {
            //  Remove possible NaN in res
            dg::blas1::copy( 0., res);
            dg::blas2::symv( m_prolongation, m_tmp, res);
        }
        else
        {
            res = m_tmp;
        }
    }
 
  private:
    std::array<unsigned, 1> coo2axis( enum coo2d direction)
    {
        if( direction == coo2d::x)
            return {0};
        return {1};
    }
    std::array<unsigned, 1> coo2axis( enum coo3d direction)
    {
        if( direction == coo3d::x)
            return {0};
        if( direction == coo3d::y)
            return {1};
        return {2};
    }
    ContainerType m_tmp;
    IMatrix m_average, m_prolongation;
    double m_area_inv;
};
 
}//namespace dg