     std::cout << "relative difference in volume is "<<fabs(volumeRZP - volume)/volume<<std::endl;
     std::cout << "Note that the error might also come from the volume in RZP!\n";
 
-    /////////////////////////TEST 3d grid//////////////////////////////////////
-    std::cout << "Start DS test!"<<std::endl;
-    const dg::DVec vol3d = dg::create::volume( g3d);
-    DFA fieldaligned( orthogonal::XField( gp, g2d, g2d.g()), g3d, gp.rk4eps, dg::NoLimiter(), dg::NEU); 
-
-    dg::DS<DFA, dg::Composite<dg::DMatrix>, dg::DVec> ds( fieldaligned, orthogonal::XField(gp, g2d, g2d.g()), dg::normed, dg::centered, false);
-    dg::DVec B = dg::pullback( solovev::InvB(gp), g3d), divB(B);
-    dg::DVec lnB = dg::pullback( solovev::LnB(gp), g3d), gradB(B);
-    dg::DVec gradLnB = dg::pullback( solovev::GradLnB(gp), g3d);
-    dg::blas1::pointwiseDivide( ones, B, B);
-
-    ds.centeredT( B, divB);
-    std::cout << "Divergence of B is "<<sqrt( dg::blas2::dot( divB, vol3d, divB))<<"\n";
-    ds.centered( lnB, gradB);
-    dg::blas1::axpby( 1., gradB, -1., gradLnB, gradLnB);
-    //test if topological shift was correct!!
-    //dg::blas1::pointwiseDot(cutter, gradLnB, gradLnB);
-    double norm = sqrt( dg::blas2::dot( gradB, vol3d, gradB) );
-    std::cout << "rel. error of lnB is    "<<sqrt( dg::blas2::dot( gradLnB, vol3d, gradLnB))/norm<<" (doesn't fullfill boundary conditions so it was cut before separatrix)\n";
-
-    const dg::DVec function = dg::pullback(solovev::FuncNeu(gp), g3d);
-    dg::DVec temp(function);
-    const dg::DVec derivative = dg::pullback(solovev::DeriNeu(gp), g3d);
-    ds( function, temp);
-    dg::blas1::axpby( 1., temp, -1., derivative, temp);
-    norm = sqrt( dg::blas2::dot( derivative, vol3d, derivative) );
-    std::cout << "rel. error of DS  is    "<<sqrt( dg::blas2::dot( temp, vol3d, temp))/norm<<"\n";
-    X = gradB;
-    //err = nc_put_var_double( ncid, divBID, periodify(X, g2d_periodic).data());
-    err = nc_put_var_double( ncid, divBID, X.data());
+    ///////////////////////////TEST 3d grid//////////////////////////////////////
+    //std::cout << "Start DS test!"<<std::endl;
+    //const dg::DVec vol3d = dg::create::volume( g3d);
+    //DFA fieldaligned( orthogonal::XField( gp, g2d, g2d.g()), g3d, gp.rk4eps, dg::NoLimiter(), dg::NEU); 
+
+    //dg::DS<DFA, dg::Composite<dg::DMatrix>, dg::DVec> ds( fieldaligned, orthogonal::XField(gp, g2d, g2d.g()), dg::normed, dg::centered, false);
+    //dg::DVec B = dg::pullback( solovev::InvB(gp), g3d), divB(B);
+    //dg::DVec lnB = dg::pullback( solovev::LnB(gp), g3d), gradB(B);
+    //dg::DVec gradLnB = dg::pullback( solovev::GradLnB(gp), g3d);
+    //dg::blas1::pointwiseDivide( ones, B, B);
+
+    //ds.centeredT( B, divB);
+    //std::cout << "Divergence of B is "<<sqrt( dg::blas2::dot( divB, vol3d, divB))<<"\n";
+    //ds.centered( lnB, gradB);
+    //dg::blas1::axpby( 1., gradB, -1., gradLnB, gradLnB);
+    ////test if topological shift was correct!!
+    ////dg::blas1::pointwiseDot(cutter, gradLnB, gradLnB);
+    //double norm = sqrt( dg::blas2::dot( gradB, vol3d, gradB) );
+    //std::cout << "rel. error of lnB is    "<<sqrt( dg::blas2::dot( gradLnB, vol3d, gradLnB))/norm<<" (doesn't fullfill boundary conditions so it was cut before separatrix)\n";
+
+    //const dg::DVec function = dg::pullback(solovev::FuncNeu(gp), g3d);
+    //dg::DVec temp(function);
+    //const dg::DVec derivative = dg::pullback(solovev::DeriNeu(gp), g3d);
+    //ds( function, temp);
+    //dg::blas1::axpby( 1., temp, -1., derivative, temp);
+    //norm = sqrt( dg::blas2::dot( derivative, vol3d, derivative) );
+    //std::cout << "rel. error of DS  is    "<<sqrt( dg::blas2::dot( temp, vol3d, temp))/norm<<"\n";
+    //X = gradB;
+    //err = nc_put_var_double( ncid, divBID, periodify(X, g3d_periodic).data());
+    ////err = nc_put_var_double( ncid, divBID, X.data());
     err = nc_close( ncid);
 
 

double cosine( double x) {return cos(x);}
typedef dg::FieldAligned< orthogonal::GridX3d<dg::DVec> , dg::IDMatrix, dg::DVec> DFA;

thrust::host_vector<double> periodify( const thrust::host_vector<double>& in, const dg::GridX2d& g)
{
    thrust::host_vector<double> out(g.size());
    for( unsigned i=0; i<g.Ny()-1; i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=0; j<g.Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[((i*g.n() + k)*g.Nx() + j)*g.n()+l];
    for( unsigned i=g.Ny()-1; i<g.Ny(); i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=0; j<g.inner_Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[((0*g.n() + k)*g.Nx() + j)*g.n()+l];
    //////////////////////////////////
    for( unsigned i=g.Ny()-1; i<g.Ny(); i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=g.inner_Nx(); j<g.Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[(((i-1)*g.n() + k)*g.Nx() + j)*g.n()+l];
    return out;
}


    dg::DVec temp0( g2d.size()), temp1(temp0);
    dg::DVec w3d = dg::create::weights( g2d);

    err = nc_put_var_double( ncid, coordsID[0], periodify(X, g3d_periodic).data());
    err = nc_put_var_double( ncid, coordsID[1], periodify(Y, g3d_periodic).data());
    //err = nc_put_var_double( ncid, coordsID[0], X.data());
    //err = nc_put_var_double( ncid, coordsID[1], Y.data());
    //err = nc_put_var_double( ncid, coord1D[0], g3d.rx0().data());
    //err = nc_put_var_double( ncid, coord1D[1], g3d.zx0().data());
    //err = nc_put_var_double( ncid, coord1D[2], g3d.rx1().data());
    //err = nc_put_var_double( ncid, coord1D[3], g3d.zx1().data());
    err = nc_put_var_double( ncid, coord1D[4], periodify(g3d.f_x(), g3d_periodic).data());
    //err = nc_put_var_double( ncid, coord1D[4], g3d.f_x().data());
    //err = nc_put_var_double( ncid, coordsID[2], g.z().data());

    dg::blas1::pointwiseDivide( g2d.g_yy(), g2d.g_xx(), temp0);
    dg::blas1::axpby( 1., ones, -1., temp0, temp0);
    X=temp0;
    err = nc_put_var_double( ncid, defID, periodify(X, g3d_periodic).data());
    //err = nc_put_var_double( ncid, defID, X.data());
    X = g2d.vol();
    err = nc_put_var_double( ncid, volID, periodify(X, g3d_periodic).data());
    //err = nc_put_var_double( ncid, volID, X.data());

    std::cout << "Construction successful!\n";

