 template<class Geometry, class IMatrix, class Matrix, class container>
 void Asela<Geometry, IMatrix, Matrix, container>::initializene( const container& src, container& target)
 { 
-    std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaN, target,src, p.eps_gamma);  //=ne-1 = Gamma (ni-1)  
-    if(  number[0] == invert_invgamma.get_max())
-      throw dg::Fail( p.eps_gamma);
+    if (p.tau[1] == 0.) {
+        dg::blas1::axpby( 1.,src, 0., target); //  ne-1 = N_i -1
+    } 
+    else {
+        std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaN, target,src, p.eps_gamma);  //=ne-1 = Gamma (ni-1)  
+        if(  number[0] == invert_invgamma.get_max())
+        throw dg::Fail( p.eps_gamma);
+    }
 }
 
 template<class G, class IMatrix, class M, class V>

        multi_pol[u].set_chi( multi_chi[u]);
    }

    
    old_gammaN.extrapolate( chi);
    std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaN, chi, y[1], p.eps_gamma); //chi= Gamma (Ni-1)
    old_gammaN.update( chi);
    if(  number[0] == invert_invgamma.get_max())
        throw dg::Fail( p.eps_gamma);
    //rhs
    dg::blas1::axpby( -1., y[0], 1.,chi,chi);        //chi=  Gamma (n_i-1) - (n_e-1) = Gamma n_i - n_e
    //polarisation
    old_phi.extrapolate( phi[0]);
    number = multigrid.direct_solve( multi_pol, phi[0], chi, p.eps_pol);
    old_phi.update( phi[0]);
    if(  number[0] == invert_pol.get_max())
        throw dg::Fail( p.eps_pol);   

    if (p.beta!=0.) apar[0] = induct(y); //computes a_par and needs correct npe
//     if (p.flrmode==1) invert_invgammaA(invgammaDIR,apar[1] ,apar[0] );             //chi= Gamma (Ni-1)
    if (p.flrmode==1) {
        if (p.tau[1] == 0.) {
          dg::blas1::axpby(1.0,apar[0],0.,apar[1]);
        } 
        else {
            old_gammaApar.extrapolate( apar[1]);
            std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaDIR, apar[1], apar[0] , p.eps_gamma);
            old_gammaApar.update( apar[1]);
            if(  number[0] == invert_invgamma.get_max())
                throw dg::Fail( p.eps_gamma); 
        }
    }
    if (p.flrmode==0) dg::blas1::axpby(1.0,apar[0],0.,apar[1]);

