average.h

Go to the documentation of this file.

#pragma once
 
#include "grid.h"
#include "weights.h"
#include "dg/blas1.h"
#include "dg/backend/average_dispatch.h"
#include "dg/backend/view.h"
 
namespace dg{
template<class container>
void simple_average( unsigned nx, unsigned ny, const container& in0, const container& in1, container& out)
{
    const double* in0_ptr = thrust::raw_pointer_cast( in0.data());
    const double* in1_ptr = thrust::raw_pointer_cast( in1.data());
          double* out_ptr = thrust::raw_pointer_cast( out.data());
    dg::View<const container> in0_view( in0_ptr, nx), in1_view( in1_ptr, nx);
    dg::View<container> out_view( out_ptr, nx);
    dg::blas1::pointwiseDot( in0_view, in1_view, out_view);
    for( unsigned i=1; i<ny; i++)
    {
        in0_view.construct( in0_ptr+i*nx, nx);
        in1_view.construct( in1_ptr+i*nx, nx);
        dg::blas1::pointwiseDot( 1., in0_view, in1_view, 1, out_view);
    }
}
 
template< class ContainerType>
struct Average
{
    using container_type = ContainerType;
    Average( const aTopology2d& g, enum coo2d direction, std::string mode = "exact") : m_mode(mode)
    {
        m_nx = g.Nx()*g.nx(), m_ny = g.Ny()*g.ny();
        m_w=dg::construct<ContainerType>(dg::create::weights(g, direction));
        m_temp = m_w;
        m_transpose = false;
        unsigned size1d = 0;
        if( direction == coo2d::x)
        {
            dg::blas1::scal( m_temp, 1./g.lx());
            dg::blas1::scal( m_w, 1./g.lx());
            size1d = m_ny;
            if( "simple" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
        }
        else
        {
            m_transpose = true;
            dg::blas1::scal( m_temp, 1./g.ly());
            dg::blas1::scal( m_w, 1./g.ly());
            if( "exact" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
            size1d = m_nx;
        }
        thrust::host_vector<double> t1d( size1d);
        m_temp1d = dg::construct<ContainerType>( t1d);
        if( !("exact"==mode || "simple" == mode))
            throw dg::Error( dg::Message( _ping_) << "Mode must either be exact or simple!");
 
    }
 
    Average( const aTopology3d& g, enum coo3d direction, std::string mode = "exact"): m_mode(mode)
    {
        m_w = dg::construct<ContainerType>(dg::create::weights(g, direction));
        m_temp = m_w;
        m_transpose = false;
        unsigned nx = g.nx()*g.Nx(), ny = g.ny()*g.Ny(), nz = g.nz()*g.Nz();
        if( direction == coo3d::x) {
            dg::blas1::scal( m_temp, 1./g.lx());
            dg::blas1::scal( m_w, 1./g.lx());
            m_nx = nx, m_ny = ny*nz;
            if( "simple" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
        }
        else if( direction == coo3d::z) {
            m_transpose = true;
            dg::blas1::scal( m_temp, 1./g.lz());
            dg::blas1::scal( m_w, 1./g.lz());
            m_nx = nx*ny, m_ny = nz;
            if( "exact" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
        }
        else if( direction == coo3d::xy) {
            dg::blas1::scal( m_temp, 1./g.lx()/g.ly());
            dg::blas1::scal( m_w, 1./g.lx()/g.ly());
            m_nx = nx*ny, m_ny = nz;
            if( "simple" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
        }
        else if( direction == coo3d::yz) {
            m_transpose = true;
            dg::blas1::scal( m_temp, 1./g.ly()/g.lz());
            dg::blas1::scal( m_w, 1./g.ly()/g.lz());
            m_nx = nx, m_ny = ny*nz;
            if( "exact" == mode)
                dg::transpose( m_nx, m_ny, m_temp, m_w);
        }
        else
            std::cerr << "Warning: this direction is not implemented\n";
        if(!m_transpose)
            m_temp1d = dg::construct<ContainerType>(
                thrust::host_vector<double>( m_ny,0.));
        else
            m_temp1d = dg::construct<ContainerType>(
                thrust::host_vector<double>( m_nx,0.));
        if( !("exact"==mode || "simple" == mode))
            throw dg::Error( dg::Message( _ping_) << "Mode must either be exact or simple!");
    }
    void operator() (const ContainerType& src, ContainerType& res, bool extend = true)
    {
        if( !m_transpose)
        {
            //temp1d has size m_ny
            if( "exact" == m_mode)
                dg::average( m_nx, m_ny, src, m_w, m_temp1d);
            else
            {
                dg::transpose( m_nx, m_ny, src, m_temp);
                dg::simple_average( m_ny, m_nx, m_temp, m_w, m_temp1d);
            }
            if( extend )
                dg::extend_column( m_nx, m_ny, m_temp1d, res);
            else
                res = m_temp1d;
        }
        else
        {
            //temp1d has size m_nx
            if( "exact" == m_mode)
            {
                dg::transpose( m_nx, m_ny, src, m_temp);
                dg::average( m_ny, m_nx, m_temp, m_w, m_temp1d);
            }
            else
                dg::simple_average( m_nx, m_ny, src, m_w, m_temp1d);
            if( extend )
                dg::extend_line( m_nx, m_ny, m_temp1d, res);
            else
                res = m_temp1d;
        }
    }
 
  private:
    unsigned m_nx, m_ny;
    ContainerType m_w, m_temp, m_temp1d;
    bool m_transpose;
    std::string m_mode;
};
 
 
}//namespace dg