     std::cout << "volumeXYP is "<< volume<<std::endl;
     std::cout << "volumeRZP is "<< volumeRZP<<std::endl;
     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";
 
 
     return 0;

    err = nc_create( "test.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
    int dim3d[3];
    err = file::define_dimensions(  ncid, dim3d, g.grid());
    int coordsID[2], onesID;
    err = nc_def_var( ncid, "r_xy", NC_DOUBLE, 3, dim3d, &coordsID[0]);
    err = nc_def_var( ncid, "z_xy", NC_DOUBLE, 3, dim3d, &coordsID[1]);
    err = nc_def_var( ncid, "psi", NC_DOUBLE, 3, dim3d, &onesID);

    thrust::host_vector<double> ones = dg::pullback( psip, g);
    g.grid().display();
    err = nc_put_var_double( ncid, onesID, ones.data());
    err = nc_put_var_double( ncid, coordsID[0], g.r().data());
    err = nc_put_var_double( ncid, coordsID[1], g.z().data());
    err = nc_close( ncid);

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

    //compute error in volume element
    dg::HVec temp0( g.grid().size()), temp1(temp0);
    dg::HVec w3d = dg::create::weights( g.grid());
    dg::blas1::pointwiseDot( g.g_xx(), g.g_yy(), temp0);

    dg::blas1::pointwiseDot( g.g_xx(), g.g_yy(), temp0);
    dg::blas1::pointwiseDot( g.g_xy(), g.g_xy(), temp1);
    dg::blas1::axpby( 1., temp0, -1., temp1, temp0);
    dg::blas1::pointwiseDot( temp0, g.g_pp(), temp0);
    dg::blas1::transform( temp0, temp0, dg::SQRT<double>());
    dg::blas1::pointwiseDivide( ones, temp0, temp0);
    dg::blas1::axpby( 1., temp0, -1., g.vol(), temp0);
    std::cout << "Rel Consistency  of volume is "<<sqrt(dg::blas2::dot( temp0, w3d, temp0)/dg::blas2::dot( g.vol(), w3d, g.vol()))<<"\n";

    dg::blas1::pointwiseDivide( g.r(), g.g_xx(), temp0);
