     }
     /// --------- More flux labels --------------------------------
     dg::Grid1d grid1d(psipmin<psipmax ? psipmin : psipmax,
-            psipmin<psipmax ? psipmax : psipmin, npsi ,Npsi,dg::NEU);
+            psipmin<psipmax ? psipmax : psipmin, npsi ,Npsi,dg::DIR_NEU); //inner value is always zero
     if( !gp.isToroidal())
     {
         dg::HVec rho = dg::evaluate( dg::cooX1d, grid1d);

        dg::SparseTensor<dg::HVec> metricX = gX2d.metric();
        dg::HVec volX2d = dg::tensor::volume2d( metricX);
        dg::blas1::pointwiseDot( coordsX[0], volX2d, volX2d); //R\sqrt{g}
        dg::HVec dvdzeta;
        avg_eta( volX2d, dvdzeta, false);
        dg::blas1::scal( dvdzeta, 4.*M_PI*M_PI);
        dg::Grid1d gX1d( gX2d.x0(), gX2d.x1(), npsi, Npsi, dg::NEU);
        dg::HVec X_psi_vol = dg::integrate( dvdzeta, gX1d);
        map1d.emplace_back( "dvdzeta", dvdzeta,
            "dvdzeta on X-point grid");
        map1d.emplace_back( "psi_vol", X_psi_vol,
            "Flux volume on X-point grid");


    std::vector<std::tuple<std::string, dg::HVec, std::string> > map1d;
    ///////////TEST CURVILINEAR GRID TO COMPUTE FSA QUANTITIES
    unsigned npsi = 3, Npsi = p.Npsi;//set number of psivalues (NPsi % 8 == 0)
    //Generate list of functions to evaluate
    std::vector< std::tuple<std::string, std::string, dg::geo::CylindricalFunctor >> map{
        {"Psip", "Flux function", mag.psip()},
