         else if( "upwind-advection" == m_advection)
         {
             dg::blas1::copy( 0., yp);
-            // v_x dx n
+            //  - v_x dx n
             dg::blas2::symv( -1., m_centered[1], m_psi, 0., m_v); //v_x
             dg::blas2::symv( m_forward[0], y, m_temp[1]);
             dg::blas2::symv( m_backward[0], y, m_temp[2]);
-            dg::blas1::subroutine( shu::UpwindAdvection(), yp, m_temp[1], m_temp[2], m_v);
-            // v_y dy n
+            dg::blas1::evaluate( yp, dg::minus_equals(), dg::UpwindProduct(), m_v, m_temp[2], m_temp[1]);
+            //  - v_y dy n
             dg::blas2::symv( 1., m_centered[0], m_psi, 0., m_v); //v_y
             dg::blas2::symv( m_forward[1], y, m_temp[1]);
             dg::blas2::symv( m_backward[1], y, m_temp[2]);
-            dg::blas1::subroutine( shu::UpwindAdvection(), yp, m_temp[1], m_temp[2], m_v);
+            dg::blas1::evaluate( yp, dg::minus_equals(), dg::UpwindProduct(), m_v, m_temp[2], m_temp[1]);
         }
         else if( "centered-advection" == m_advection)
         {

namespace shu
{

struct Upwind{
    DG_DEVICE
    void operator()( double& result, double fw, double bw, double v)
    {
        if( v < 0)
            result -= fw; // yp = - Div( F)
        else
            result -= bw;
    }
};
struct UpwindAdvection{
    DG_DEVICE
    void operator()( double& result, double fw, double bw, double v)
    {
        if( v < 0)
            result -= v*fw; // yp = - v Grad f
        else
            result -= v*bw;
    }
};
template< class Geometry, class Matrix, class Container>
struct Diffusion
{

            dg::blas2::symv( m_inter, m_v, m_fine_v);
            dg::blas2::symv( m_inter, m_temp[1], m_fine_temp[1]);
            dg::blas2::symv( m_inter, m_temp[2], m_fine_temp[2]);
            dg::blas1::evaluate( m_fine_yp, dg::minus_equals(), dg::UpwindProduct(), m_fine_v, m_fine_temp[2], m_fine_temp[1]);
            // v_y dy n
            dg::blas2::symv( 1., m_centered[0], m_psi, 0., m_v); //v_y
            dg::blas2::symv( m_forward[1], y, m_temp[1]);
