             for (unsigned zz=0;zz<1;zz++) //Nz iterator
             {
                 std::cout << "n = " << k*n << " Nx = " <<pow(2,i)* Nx << " Ny = " <<pow(2,i)* Ny << " Nz = "<<pow(2,zz)* Nz <<"\n";
-                dg::Grid3d<double> g3d( Rmin,Rmax, Zmin,Zmax, 0, 2.*M_PI,k*n,pow(2,i)* Nx,pow(2,i)* Ny, pow(2,zz)*Nz);
+                dg::Grid3d<double> g3d( Rmin,Rmax, Zmin,Zmax, 0, 2.*M_PI,k*n,pow(2,i)* Nx,pow(2,i)* Ny, pow(2,zz)*Nz,dg::NEU, dg::NEU, dg::PER, dg::cylindrical);
                 const dg::DVec w3d = dg::create::weights( g3d);
                 dg::DVec pupilongrid = dg::evaluate( pupil, g3d);
 

    err = nc_put_var_double( ncid, vecID[1], vecZ.data());
    err = nc_put_var_double( ncid, vecID[2], vecP.data());
    nc_close(ncid);
    
    dg::HVec v5(1, 0);
    std::vector<thrust::host_vector<double> > in(3, v5);
    std::vector<thrust::host_vector<double> > out(3, v5);
    in[0][0]=gp.R_0+0.9*gp.a; 
    in[1][0]=0.9*gp.a*gp.elongation;
    in[2][0]=0.;

    dg::integrateRK4( field, in, out,  2*M_PI, gp.rk4eps);
    
    std::cout << "Check single field by integrating from 0 to 2pi" << "\n";
    std::cout <<"Rin =  "<< in[0][0] <<" Zin =  "<<in[1][0] <<" sin  = "<<in[2][0]<<"\n";
    std::cout <<"Rout = "<< out[0][0]<<" Zout = "<<out[1][0]<<" sout = "<<out[2][0]<<"\n";



    std::vector<double> c(13);
    //read and store geom data
       std::vector<double> v, v2;
    try{ 
        if( argc==1)
        {
            v = file::read_input( "geometry_params.txt"); 
            v2 = file::read_input( "../feltor/input.txt");

        }
        else
        {
            v = file::read_input( argv[1]); 
            v2 = file::read_input( argv[2]);
        }
    }
    catch (toefl::Message& m) {  
        m.display(); 
        for( unsigned i = 0; i<v.size(); i++)
            std::cout << v[i] << " ";
            std::cout << std::endl;
        return -1;}

    //write parameters from file into variables
    const solovev::GeomParameters gp(v);
    const eule::Parameters p(v2);
    p.display( std::cout);
    gp.display( std::cout);

    Rmin=gp.R_0-(p.boxscale)*gp.a;
    Zmin=-(p.boxscale)*gp.a*gp.elongation;
    Rmax=gp.R_0+(p.boxscale)*gp.a; 
    Zmax=(p.boxscale)*gp.a*gp.elongation;
    std::cout << "The grid parameters" <<"\n";
    std::cout  << Rmin<<" rho_s " << Rmax <<" rho_s " << Zmin <<" rho_s " <<Zmax <<" rho_s " <<"\n";
    std::cout << "Type n, Nx, Ny, Nz\n";
    std::cin >> n>> Nx>>Ny>>Nz;
        
    solovev::Field field(gp);
    solovev::Psip psip(gp.R_0,gp.A,gp.c);
    solovev::PsipR psipR(gp.R_0,gp.A,gp.c);
