 
 typedef cusp::device_memory Memory;
 
+double initial( double x, double y){ return 2.*sin(x)*sin(y);}
+double solution( double x, double y){ return 2.*sin(x)*sin(y)*exp(-2.*D*T);}
 using namespace std;
 
 int main()

    cout << "# of grid cells:            " << Nx*Ny<<endl;
    cout << "Timestep                    " << dt << endl;
    cout << "Diffusion                   " << D <<endl;
    dg::Lamb lamb( 0.5*lx, 0.5*ly, R, U);
    HArrVec omega = expand< dg::Lamb, n> ( lamb, 0, lx, 0, ly, Nx, Ny);
    DArrVec stencil = expand< double(&)(double, double), n> ( one, 0, lx, 0, ly, Nx, Ny);
    dg::Lamb lamb2( 0.5*lx, 0.5*ly-0.9755*U*T, R, U);
    HArrVec solh = expand< dg::Lamb, n> ( lamb2, 0, lx, 0, ly, Nx, Ny);
    DVec sol = solh.data();
    DVec y0( omega.data()), y1( y0);
    //make solver and stepper

        time += dt;
    }
    ////////////////////////////////////////////////////////////////////
    //cout << "Analytic formula enstrophy "<<lamb.enstrophy()<<endl;
    //cout << "Analytic formula energy    "<<lamb.energy()<<endl;
    cout << "Total vorticity           is: "<<blas2::dot( stencil.data(), S2D<double, n>(hx, hy), y0) << "\n";
    cout << "Relative enstrophy error  is: "<<(0.5*blas2::dot( S2D<double, n>(hx, hy), y0) - enstrophy)/enstrophy<<"\n";
    test( y0, y1); //get the potential ready
    cout << "Relative energy error     is: "<<(0.5*blas2::dot( test.potential(), S2D<double, n>(hx, hy), y0) - energy)/energy<<"\n";

    blas1::axpby( 1., sol, -1., y0);
    cout << "Distance to solution "<<sqrt( blas2::dot( S2D<double,n >(hx,hy), y0)) << endl;

    //energy and enstrophy errrors are due to timestep only ( vorticity is exactly conserved)
    // k = 2 | p = 3
    // k = 3 | p = 4
    // k = 4 | p = 5

    //solution to sin(x)sin(y) 
    // n = 1 
    // n = 2 | p = 2
    // n = 3 | p = 2.6
    // n = 4 | p = 4

    return 0;

}
