      * @return phi[0] is the electron and phi[1] the generalized ion potential
      */
     const std::vector<container>& potential( ) const { return phi;}
-    #if def APAR
-        const container& aparallel( ) const { return apar;}
-    #endif
 
     /**
      * @brief Return the Gamma operator used by this object

    double Uphi = 0.5*p.mu[1]*dg::blas2::dot( expy[1], w3d, omega); 
    dg::blas1::pointwiseDot( y[2], y[2], omega); //U_e^2
    double Upare = -0.5*p.mu[0]*dg::blas2::dot( expy[0], w3d, omega); 
    dg::blas1::pointwiseDot( y[3], y[3], omega); //U_i^2
    double Upari =  0.5*p.mu[1]*dg::blas2::dot( expy[1], w3d, omega); 
    energy_ = Ue + Ui + Uphi + Upare + Upari;

    for( unsigned i=0; i<2; i++)
    {

        arakawa( y[i], phi[i], yp[i]);                      
        arakawa( y[i+2], phi[i], yp[i+2]);  
#if def APAR
        arakawa( apar, y[i], arakAN[i]); // [A_parallel,N]_RZ
        arakawa( apar, u[i], arakAU[i]); // [A_parallel,U]_RZ
#endif
        dg::blas1::pointwiseDot( yp[i], binv, yp[i]);                        // dtlnN =1/B [phi,lnN]_RZ
        dg::blas1::pointwiseDot( yp[2+i], binv, yp[2+i]);                    // dtU =1/B [phi,U]_RZ
    
        //compute parallel derivatives
        dz(y[i], dzy[i]);
        dz(phi[i], dzphi[i]);
        dz(y[2+i], dzy[2+i]);
    
        //add A_parallel terms to the parallel derivaitives
//         dg::axpby(-p.beta,araAparaN ,1.,dzy[i]);
//         dg::axpby(-p.beta,araAparaphi ,1.,dzphi[i]);
//         dg::axpby(-p.beta,araAparaU,1.,dzy[2+i]);
        //parallel advection terms
        dg::blas1::pointwiseDot(y[2+i], dzy[i], omega); //Udz lnN 
        dg::blas1::axpby( -1., omega, 1., yp[i]);                            //dtlnN = dtlnN - Udz lnN
        dg::blas1::axpby( -1., dzy[2+i], 1., yp[i]);                         //dtlnN = dtlnN - dz U
        dg::blas1::pointwiseDot(y[2+i], gradlnB, omega);                     //dtlnN = dtlnN + U dz ln B
        dg::blas1::axpby( 1., omega, 1., yp[i]); 
        dg::blas1::pointwiseDot(y[2+i], dzy[2+i], omega);                    //dtU = dtU - U dz U
        dg::blas1::axpby( -1., omega, 1., yp[2+i]);

        //parallel force terms
        dg::blas1::axpby( -p.tau[i]/p.mu[i]/p.eps_hat, dzy[i], 1., yp[2+i]); //dtU = dtU - tau/(hat(mu))*dz lnN
        dg::blas1::axpby( -1./p.mu[i]/p.eps_hat, dzphi[i], 1., yp[2+i]);     //dtU = dtU - 1/(hat(mu))*dz phi

        //curvature terms
        curve( y[i], curvy[i]); //K(N)
        curve( y[2+i], curvy[2+i]);//K(U)
        curve( phi[i], curvphi[i]);//K(phi)

        dg::blas1::pointwiseDot( y[2+i], curvy[2+i], omega); //U K(U)
        dg::blas1::pointwiseDot( y[2+i], omega, chi); //U^2 K(U)
        dg::blas1::axpby( -p.mu[i]*p.eps_hat, omega, 1., yp[i]);             //dtlnN = dtlnN - (hat(mu)) U K(U)
        dg::blas1::axpby( -0.5*p.mu[i]*p.eps_hat, chi, 1., yp[2+i]);         //dtU = dtU - 0.5 (hat(mu)) U^2 K(U)

        dg::blas1::pointwiseDot( y[2+i], curvy[i], omega); //U K(ln N)
        dg::blas1::pointwiseDot( y[2+i], omega, chi); //U^2K(ln N)
        dg::blas1::axpby( -p.tau[i], omega, 1., yp[2+i]);                    //dtU = dtU - tau U K(lnN)
        dg::blas1::axpby( -0.5*p.mu[i]*p.eps_hat, chi, 1., yp[i]);           //dtlnN = dtlnN -0.5mu U^2K(lnN)

        dg::blas1::axpby( -p.tau[i], curvy[i], 1., yp[i]);                   //dtlnN = dtlnN - tau K(lnN)
        dg::blas1::axpby( -2.*p.tau[i], curvy[2+i], 1., yp[2+i]);            //dtU = dtU - 2 tau K(U)
        dg::blas1::axpby( -1., curvphi[i], 1., yp[i]);                       //dtlnN= dtlnN - K(psi)

        dg::blas1::pointwiseDot( y[2+i], curvphi[i], omega);  //U K(phi)
        dg::blas1::axpby( -0.5, omega, 1., yp[2+i]);                         //dtU = dtU -0.5 UK(psi)


    }
    //add parallel resistivity
//     dg::blas1::pointwiseDot( expy[0], y[2], omega);
//     dg::blas1::pointwiseDot( expy[1], y[3], chi);
//     dg::blas1::axpby( -1., omega, 1., chi); //-N_eU_e + N_iU_i
//     dg::blas1::pointwiseDivide( chi, expy[0], omega);//J_par/N_e
//     dg::blas1::pointwiseDivide( chi, expy[1], chi); //J_par/N_i
// 
//     dg::blas1::axpby( -p.c/p.mu[0]/p.eps_hat, omega, 1., yp[2]);
//     dg::blas1::axpby( -p.c/p.mu[1]/p.eps_hat, chi, 1., yp[3]);

    //add parallel diffusion with naive implementation
    for( unsigned i=0; i<4; i++)
    {

    for( unsigned i=0; i<2; i++)
    {
        dg::blas1::pointwiseDivide( source, expy[i], omega); //source/N
        dg::blas1::axpby( 1., omega, 1, yp[i]  );       //dtlnN = dtlnN + source/N
    }

    for( unsigned i=0; i<4; i++) //damping and pupil on N and w
    {
        dg::blas1::pointwiseDot( damping, yp[i], yp[i]); 
        dg::blas1::pointwiseDot( iris, yp[i], yp[i]);
