 //     ellipticsym.symv(function,dsTds);
 //     dg::blas1::scal(dsTds,-1.0);
 // //     ds.centeredT(ones,divbT);
-    ds.forwardT( derivativef, dsTdsf);  //dsT(ds(f))
+//     ds.forwardT( derivativef, dsTdsf);  //dsT(ds(f))
 //     ds.backwardT( derivativeb, dsTdsb); //dsT(ds(f))
 
 //     //centered

    const dg::DVec w2d = dg::create::weights( g2d);
    const dg::DVec v3d = dg::create::inv_volume( g3d);

    std::cout << "computing dsDIR" << std::endl;
    dg::DDS::FieldAligned dsFA( field, g3d_fein, rk4eps, dg::DefaultLimiter(), dg::DIR);
    std::cout << "computing dsNEU" << std::endl;
    dg::DDS::FieldAligned dsNUFA( field, g3d_fein, rk4eps, dg::DefaultLimiter(), dg::NEU);

    dg::DDS ds ( dsFA, g3d, field, dg::not_normed, dg::centered), 
        dsNU ( dsNUFA, g3d, field, dg::not_normed, dg::centered);

//     dg::DS<dg::DMatrix, dg::DVec> dsNEU( field, g3d, g3d.hz(), rk4eps, dg::DefaultLimiter(), dg::NEU);
    
//     dg::Grid3d<double> g3dp( Rmin,Rmax, Zmin,Zmax, z0, z1,  n, Nx, Ny, 1);

//     
//     
  
//     dsNU( function, derivative); //ds(f)

//     dsNU.forward( function, derivativef); //ds(f)
//     dsNU.backward( function, derivativeb); //ds(f)
