     double z0 = 0, z1 = 2.*M_PI;
     std::cout << "Type n, Nx, Ny, Nz\n";    
     std::cin >> n >> Nxn >> Nyn >> Nzn;
-    
+
     dg::CylindricalGrid3d g3d( Rmin,Rmax, Zmin,Zmax, z0, z1,  n,Nxn ,Nyn, Nzn,dg::DIR, dg::DIR, dg::PER);
     dg::Grid2d g2d( Rmin,Rmax, Zmin,Zmax,  n, Nxn ,Nyn);
     

    
    std::cout << "computing dsDIR" << std::endl;
    dg::geo::Fieldaligned<dg::aProductGeometry3d, dg::IDMatrix, dg::DVec>  dsFA( mag, g3d, dg::DIR, dg::DIR, dg::geo::FullLimiter(), rk4eps, mx, my, true, true, integrate);
    std::cout << "computing dsNEU" << std::endl;
    dg::geo::Fieldaligned<dg::aProductGeometry3d, dg::IDMatrix, dg::DVec> dsNUFA( mag, g3d,dg::NEU, dg::NEU, dg::geo::FullLimiter(), rk4eps, mx, my, true, true, integrate);

    dg::geo::DS<dg::aProductGeometry3d, dg::IDMatrix, dg::DMatrix, dg::DVec> ds ( dsFA, dg::not_normed, dg::centered), 
        dsNU ( dsNUFA, dg::not_normed, dg::centered);
    
    dg::DVec function = dg::evaluate( funcNEU, g3d) ,
                        temp( function),


    const dg::DVec gradlnB_ = dg::evaluate(gradlnB, g3d);
  
    ds( function, derivative); //ds(f)
    dg::blas1::pointwiseDivide(ones,  inverseB, temp); //B
    ds.centeredDiv( temp, divBT); // dsT B
    ds.centeredDiv( ones, divbT);

    double norm = dg::blas2::dot( w3d, solution);
    double err =dg::blas2::dot( w3d, derivative);
    dg::blas1::axpby( 1., solution, -1., derivative);
    err =dg::blas2::dot( w3d, derivative);
    std::cout << "Relative Difference in ds f   = "<< sqrt( err/norm )<<"\n"; 

    norm = dg::blas2::dot( w3d, divbsol);
    err = dg::blas2::dot( w3d, divbT);
    dg::blas1::axpby( 1., divbsol, -1., divbT);
    err = dg::blas2::dot(divbT, w3d,divbT);
    std::cout << "Relative Difference in div(b) = "<<sqrt(err/norm)<<"\n";
    double normdivBT =dg::blas2::dot(divBT, w3d,divBT);
    std::cout << "Error in div(B)  = "<<sqrt( normdivBT)<<"\n";  
    
    return 0;
}
