     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);
+    try{ v3 = file::read_input( "geometry_params.txt"); }
+    catch (toefl::Message& m) {  
+        m.display(); 
+        for( unsigned i = 0; i<v.size(); i++)
+//             std::cout << v3[i] << " ";
+//             std::cout << std::endl;
+        return -1;}
+
+    const solovev::GeomParameters gp(v3);
+    gp.display( std::cout);
+    double Rmin=gp.R_0-(1.05)*gp.a;
+    double Zmin=-(1.05)*gp.a*gp.elongation;
+    double Rmax=gp.R_0+(1.05)*gp.a; 
+    double Zmax=(1.05)*gp.a*gp.elongation;
+    //Make grid
+     dg::Grid3d<double > grid( Rmin,Rmax, Zmin,Zmax, 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);
+    eule::Feltor< dg::DVec > feltor( grid, p,gp); 
+    eule::Rolkar< dg::DVec > rolkar( grid, p,gp);
 
     //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);
+    //with bath
+      dg::BathRZ init0(16,16,p.Nz,Rmin,Zmin, 30.,15.,p.amp);
+    //with gaussians
+//     dg::Gaussian3d init0( p.R_0, p.posY*p.a,    M_PI, p.sigma, p.sigma, M_PI/8.*p.m_par, p.amp );     
+//     dg::Gaussian3d init1( p.R_0, -p.a*p.posY,   M_PI, p.sigma, p.sigma, M_PI/8.*p.m_par, p.amp ); 
+//     dg::Gaussian3d init2( p.R_0+p.posX*p.a, 0., M_PI, p.sigma, p.sigma, M_PI/8.*p.m_par, p.amp ); 
+//     dg::Gaussian3d init3( p.R_0-p.a*p.posX, 0., M_PI, p.sigma, p.sigma, M_PI/8.*p.m_par, p.amp ); 
+//     eule::Gradient grad(p.R_0, p.a, p.a-p.b, p.lnn_inner);
+    solovev::Gradient grad(gp, p.lnn_inner);
+  
+    //const dg::HVec gradient( dg::evaluate(grad, grid));
+    std::vector<dg::DVec> y0(4, dg::evaluate( grad, grid)), y1(y0); 
+    dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init0, 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., (dg::DVec)dg::evaluate(init0, grid), 1., y0[1]);
+//     dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init1, grid), 1., y0[1]);
+//     dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init2, grid), 1., y0[1]);
+//     dg::blas1::axpby( 1., (dg::DVec)dg::evaluate(init3, grid), 1., y0[1]);
 
-    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::blas1::pointwiseDot(rolkar.iris(),y0[0],y0[0]);
+    dg::blas1::pointwiseDot(rolkar.iris(),y0[1],y0[1]);
     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.h"
#include "dg/timer.cuh"



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