evaluation.h

Go to the documentation of this file.

#pragma once
 
#include <cassert>
#include <cmath>
#include <thrust/host_vector.h>
#include "dg/backend/config.h"
#include "grid.h"
#include "operator.h"
 
namespace dg
{
 
template< class Functor, class Topology, size_t ...I>
auto do_evaluate( Functor&& f, const Topology& g, std::index_sequence<I...>)
{
    return dg::kronecker( std::forward<Functor>(f), g.abscissas(I)...);
}
 
 
template< class Functor, class Topology>
auto evaluate( Functor&& f, const Topology& g)
{
    // The evaluate function is the reason why our Topology needs to have fixed
    // sized dimensions instead of dynamically sized dimensions
    // even though if we really wanted we could maybe ask ndim = g.ndim()
    // and then do use switch and manually implement until ndim < 10 say
    // for now we keep fixed sized grids ...
    //
    // If we ever change the order of the fastest dimension we need to rethink
    // NetCDF hyperslab construction
    return do_evaluate( std::forward<Functor>(f), g, std::make_index_sequence<Topology::ndim()>());
};
 
//These overloads help the compiler in a situation where a free function has
//several overloads of different dimensions e.g.
//double zero( double);
//double zero( double,double);
//In such a case dg::evaluate( zero, grid); cannot determine the Functor type...
template<class Topology, class value_type = typename Topology::value_type, class result_type = typename Topology::value_type, typename = std::enable_if_t<Topology::ndim() == 1 > >
auto evaluate( result_type (*f)( value_type), const Topology& g)
{
    return do_evaluate( f, g, std::make_index_sequence<Topology::ndim()>());
}
template<class Topology, class value_type0 = typename Topology::value_type, class value_type1 = typename Topology::value_type, class result_type = typename Topology::value_type, typename = std::enable_if_t<Topology::ndim() == 2 > >
auto evaluate( result_type (*f)( value_type0, value_type1), const Topology& g)
{
    return do_evaluate( f, g, std::make_index_sequence<Topology::ndim()>());
}
template<class Topology, class value_type0 = typename Topology::value_type, class value_type1 = typename Topology::value_type, class value_type2 = typename Topology::value_type, class result_type = typename Topology::value_type, typename = std::enable_if_t<Topology::ndim() == 3 > >
auto evaluate( result_type (*f)( value_type0, value_type1, value_type2), const Topology& g)
{
    return do_evaluate( f, g, std::make_index_sequence<Topology::ndim()>());
}
 
 
 
 
template<class real_type>
thrust::host_vector<real_type> integrate( const thrust::host_vector<real_type>& in, const RealGrid<real_type,1>& g, dg::direction dir = dg::forward)
{
    double h = g.hx();
    unsigned n = g.nx();
    thrust::host_vector<real_type> to_out(g.size(), 0.);
    thrust::host_vector<real_type> to_in(in);
    if( dir == dg::backward ) //reverse input vector
    {
        for( unsigned i=0; i<in.size(); i++)
            to_in[i] = in[ in.size()-1-i];
    }
 
 
    dg::SquareMatrix<real_type> forward = dg::DLT<real_type>::forward(n);
    dg::SquareMatrix<real_type> backward = dg::DLT<real_type>::backward(n);
    dg::SquareMatrix<real_type> ninj = create::ninj<real_type>( n );
    SquareMatrix<real_type> t = create::pipj_inv<real_type>(n);
    t *= h/2.;
    ninj = backward*t*ninj*forward;
    real_type constant = 0.;
 
    for( unsigned i=0; i<g.Nx(); i++)
    {
        for( unsigned k=0; k<n; k++)
        {
            for( unsigned l=0; l<n; l++)
                to_out[ i*n + k] += ninj(k,l)*to_in[ i*n + l];
            to_out[ i*n + k] += constant;
        }
        for( unsigned l=0; l<n; l++)
            constant += h*forward(0,l)*to_in[i*n+l];
    }
    thrust::host_vector<real_type> out(to_out);
    if( dir == dg::backward ) //reverse output
    {
        for( unsigned i=0; i<in.size(); i++)
            out[i] = -to_out[ in.size()-1-i]; // minus from reversing!
    }
    return out;
}
 
 
template< class UnaryOp,class real_type>
thrust::host_vector<real_type> integrate( UnaryOp f, const RealGrid<real_type,1>& g, dg::direction dir = dg::forward)
{
    thrust::host_vector<real_type> vector = evaluate( f, g);
    return integrate<real_type>(vector, g, dir);
}
template<class real_type>
thrust::host_vector<real_type> integrate( real_type (f)(real_type), const RealGrid<real_type,1>& g, dg::direction dir = dg::forward)
{
    thrust::host_vector<real_type> vector = evaluate( f, g);
    return integrate<real_type>(vector, g, dir);
};
 
}//namespace dg