 
 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()

using namespace std;
using namespace dg;

const unsigned n = 3;
const unsigned Nx = 33;
const unsigned Ny = 33;
const double lx = 1.;
const double ly = 1.;

const unsigned k = 3;
const double D = 0.0;

    cout << "# of grid cells:            " << Nx*Ny<<endl;
    cout << "Timestep                    " << dt << endl;
    cout << "Diffusion                   " << D <<endl;
    DArrVec stencil = expand< double(&)(double, double), n> ( one, 0, lx, 0, ly, Nx, Ny);

    //dg::Lamb lamb( 0.5*lx, 0.5*ly, R, U);
    //HArrVec omega = expand< dg::Lamb, n> ( lamb, 0, lx, 0, ly, Nx, Ny);
    HArrVec omega = expand< double(&)(double, double), n> ( initial, 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);
    HArrVec solh = expand< double(&)(double, double), n> ( solution, 0, lx, 0, ly, Nx, Ny);

    DVec sol = solh.data();
    DVec y0( omega.data()), y1( y0);
    //make solver and stepper
