 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::blas1::pointwiseDot( m_fine_y, m_fine_v, m_fine_temp[0]); //f_x
            dg::blas2::symv( m_fine_forward[0], m_fine_temp[0], m_fine_temp[1]);
            dg::blas2::symv( m_fine_backward[0], m_fine_temp[0], m_fine_temp[2]);
            dg::blas1::subroutine( shu::Upwind(), m_fine_yp, m_fine_temp[1], m_fine_temp[2], m_fine_v);
            //dy ( nv_y)
            dg::blas2::symv( 1., m_fine_centered[0], m_fine_psi, 0., m_fine_v); //v_y
            dg::blas1::pointwiseDot( m_fine_y, m_fine_v, m_fine_temp[0]); //f_y
            dg::blas2::symv( m_fine_forward[1], m_fine_temp[0], m_fine_temp[1]);
            dg::blas2::symv( m_fine_backward[1], m_fine_temp[0], m_fine_temp[2]);
            dg::blas1::subroutine( shu::Upwind(), m_fine_yp, m_fine_temp[1], m_fine_temp[2], m_fine_v);
        }
        else if( "centered" == m_advection)
        {

        else if( "upwind" == m_advection)
        {
            dg::blas1::copy( 0., yp);
            // - dx ( nv_x)
            dg::blas2::symv( -1., m_centered[1], m_psi, 0., m_v); //v_x
            dg::blas1::pointwiseDot( y, m_v, m_temp[0]); //f_x
            dg::blas2::symv( m_forward[0], m_temp[0], m_temp[1]);
            dg::blas2::symv( m_backward[0], m_temp[0], m_temp[2]);
            dg::blas1::evaluate( yp, dg::minus_equals(), dg::Upwind(), m_v, m_temp[2], m_temp[1]);
            // - dy ( nv_y)
            dg::blas2::symv( 1., m_centered[0], m_psi, 0., m_v); //v_y
            dg::blas1::pointwiseDot( y, m_v, m_temp[0]); //f_y
            dg::blas2::symv( m_forward[1], m_temp[0], m_temp[1]);
            dg::blas2::symv( m_backward[1], m_temp[0], m_temp[2]);
            dg::blas1::evaluate( yp, dg::minus_equals(), dg::Upwind(), m_v, m_temp[2], m_temp[1]);
        }
        else if( "centered" == m_advection)
        {
            // - dx ( nv_x)
            dg::blas2::symv( -1., m_centered[1], m_psi, 0., m_v); //v_x
            dg::blas1::pointwiseDot( y, m_v, m_temp[0]); //f_x
            dg::blas2::symv( -1., m_centered[0], m_temp[0], 0., yp);
            // - dy ( nv_y)
            dg::blas2::symv( 1., m_centered[0], m_psi, 0., m_v); //v_y
            dg::blas1::pointwiseDot( y, m_v, m_temp[0]); //f_y
            dg::blas2::symv( -1., m_centered[1], m_temp[0], 1., yp);
