 int main( int argc, char* argv[])
 {
     //Parameter initialisation
-    std::vector<double> v;
+    std::vector<double> v,v3;
     std::string input;
     if( argc != 3)
     {

    p.display( std::cout);

    ////////////////////////////////set up computations///////////////////////////
    dg::Grid3d<double > grid( p.R_0-p.a*(1+1e-1), p.R_0 + p.a*(1+1e-1),  -p.a*(1+1e-1), p.a*(1+1e-1), 0, 2.*M_PI, p.n, p.Nx, p.Ny, p.Nz, dg::DIR, dg::DIR, dg::PER);
    //create RHS 
    eule::Feltor< dg::DVec > feltor( grid, p); 
    eule::Rolkar< dg::DVec > rolkar( grid, p.nu_perp, p.nu_parallel,p.R_0, p.a, p.b, p.mu[0]*p.eps_hat);

    //create initial vector
    dg::Gaussian3d init0( p.R_0, p.posY*p.a,    0, p.sigma, p.sigma, M_PI/8., p.amp ); 
    dg::Gaussian3d init1( p.R_0, -p.a*p.posY,   0, p.sigma, p.sigma, M_PI/8., p.amp ); 
    dg::Gaussian3d init2( p.R_0+p.posX*p.a, 0., 0.,p.sigma, p.sigma, M_PI/8., p.amp ); 
    dg::Gaussian3d init3( p.R_0-p.a*p.posX, 0., 0.,p.sigma, p.sigma, M_PI/8., p.amp ); 
    eule::Gradient grad(p.R_0, p.a, p.a-p.b, p.lnn_inner);

    std::vector<dg::DVec> y0(4, dg::evaluate( init0, grid)); // n_e' = gaussian
    dg::blas1::axpby( 1., (dg::DVec)dg::evaluate( grad, grid), 1., y0[0]);
    dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init1, grid), 1., y0[0]);
    dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init2, grid), 1., y0[0]);
    dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init3, grid), 1., y0[0]);

    dg::blas1::axpby( 1., y0[0], 0., y0[1]);
    dg::blas1::axpby( 0., y0[2], 0., y0[2]); //set U = 0
    dg::blas1::axpby( 0., y0[3], 0., y0[3]); //set U = 0

    feltor.log( y0, y0, 2); //transform to logarithmic values
    dg::Karniadakis< std::vector<dg::DVec> > karniadakis( y0, y0[0].size(), p.eps_time);
    karniadakis.init( feltor, rolkar, y0, p.dt);
    double time = 0;

#include "file/read_input.h"

#include "parameters.h"
#include "geometry_circ.h"
#include "dg/timer.cuh"


