fieldaligned.h

Go to the documentation of this file.

#pragma once
#include <cmath>
#include <array>
 
#include "dg/algorithm.h"
#include "magnetic_field.h"
#include "fluxfunctions.h"
#include "curvilinear.h"
 
namespace dg{
namespace geo{
 
enum whichMatrix
{
    einsPlus = 0,   
    einsPlusT, 
    einsMinus, 
    einsMinusT,
    zeroPlus,  
    zeroMinus, 
    zeroPlusT, 
    zeroMinusT, 
    zeroForw  
};
 
typedef ONE FullLimiter;
 
typedef ZERO NoLimiter;
namespace detail{
 
template<class Container>
struct Container2EllSparseBlockMat;
 
// https://stackoverflow.com/questions/5052211/changing-value-type-of-a-given-stl-container
template<class real_type, class ... Args, template <class...> class Vector>
struct Container2EllSparseBlockMat<Vector<real_type, Args...>>
{
    using type = EllSparseBlockMat<real_type, Vector>;
};
 
static inline void parse_method( std::string method, std::string& i, std::string& p, std::string& f)
{
    f = "dg";
    if( method == "dg")                 i = "dg",       p = "dg";
    else if( method == "linear")        i = "linear",   p = "dg";
    else if( method == "cubic")         i = "cubic",    p = "dg";
    else if( method == "nearest")       i = "nearest",  p = "dg";
    else if( method == "dg-nearest")      i = "dg",      p = "nearest";
    else if( method == "linear-nearest")  i = "linear",  p = "nearest";
    else if( method == "cubic-nearest")   i = "cubic",   p = "nearest";
    else if( method == "nearest-nearest") i = "nearest",  p = "nearest";
    else if( method == "dg-linear")      i = "dg",       p = "linear";
    else if( method == "linear-linear")  i = "linear",   p = "linear";
    else if( method == "cubic-linear")   i = "cubic",    p = "linear";
    else if( method == "nearest-linear") i = "nearest",  p = "linear";
    else if( method == "dg-equi")       i = "dg",       p = "dg", f = "equi";
    else if( method == "linear-equi")   i = "linear",   p = "dg", f = "equi";
    else if( method == "cubic-equi")    i = "cubic",    p = "dg", f = "equi";
    else if( method == "nearest-equi")  i = "nearest",  p = "dg", f = "equi";
    else if( method == "dg-equi-nearest")      i = "dg",      p = "nearest", f = "equi";
    else if( method == "linear-equi-nearest")  i = "linear",  p = "nearest", f = "equi";
    else if( method == "cubic-equi-nearest")   i = "cubic",   p = "nearest", f = "equi";
    else if( method == "nearest-equi-nearest") i = "nearest", p = "nearest", f = "equi";
    else if( method == "dg-equi-linear")      i = "dg",       p = "linear", f = "equi";
    else if( method == "linear-equi-linear")  i = "linear",   p = "linear", f = "equi";
    else if( method == "cubic-equi-linear")   i = "cubic",    p = "linear", f = "equi";
    else if( method == "nearest-equi-linear") i = "nearest",  p = "linear", f = "equi";
    else
        throw Error( Message(_ping_) << "The method "<< method << " is not recognized\n");
}
 
struct DSFieldCylindrical3
{
    DSFieldCylindrical3( const dg::geo::CylindricalVectorLvl0& v): m_v(v){}
    void operator()( double, const std::array<double,3>& y,
            std::array<double,3>& yp) const {
        double R = y[0], Z = y[1];
        double vz = m_v.z()(R, Z);
        yp[0] = m_v.x()(R, Z)/vz;
        yp[1] = m_v.y()(R, Z)/vz;
        yp[2] = 1./vz;
    }
    private:
    dg::geo::CylindricalVectorLvl0 m_v;
};
 
struct DSFieldCylindrical4
{
    DSFieldCylindrical4( const dg::geo::CylindricalVectorLvl1& v): m_v(v){}
    void operator()( double, const std::array<double,3>& y,
            std::array<double,3>& yp) const {
        double R = y[0], Z = y[1];
        double vx = m_v.x()(R,Z);
        double vy = m_v.y()(R,Z);
        double vz = m_v.z()(R,Z);
        double divvvz = m_v.divvvz()(R,Z);
        yp[0] = vx/vz;
        yp[1] = vy/vz;
        yp[2] = divvvz*y[2];
    }
 
    private:
    dg::geo::CylindricalVectorLvl1 m_v;
};
 
struct DSField
{
    DSField() = default;
    //z component of v may not vanish
    DSField( const dg::geo::CylindricalVectorLvl1& v,
            const dg::aGeometry2d& g ):
        m_g(g)
    {
        dg::HVec v_zeta, v_eta;
        dg::pushForwardPerp( v.x(), v.y(), v_zeta, v_eta, g);
        dg::HVec vx = dg::pullback( v.x(), g);
        dg::HVec vy = dg::pullback( v.y(), g);
        dg::HVec vz = dg::pullback( v.z(), g);
        dg::HVec divvvz = dg::pullback( v.divvvz(), g);
        dg::blas1::pointwiseDivide(v_zeta,  vz, v_zeta);
        dg::blas1::pointwiseDivide(v_eta,   vz, v_eta);
        dzetadphi_  = dg::forward_transform( v_zeta, g );
        detadphi_   = dg::forward_transform( v_eta, g );
        dvdphi_     = dg::forward_transform( divvvz, g );
    }
    //interpolate the vectors given in the constructor on the given point
    void operator()(double, const std::array<double,3>& y, std::array<double,3>& yp) const
    {
        // shift point into domain
        yp[0] = interpolate(dg::lspace, dzetadphi_, y[0], y[1], *m_g);
        yp[1] = interpolate(dg::lspace, detadphi_,  y[0], y[1], *m_g);
        yp[2] = interpolate(dg::lspace, dvdphi_,    y[0], y[1], *m_g)*y[2];
    }
    private:
    thrust::host_vector<double> dzetadphi_, detadphi_, dvdphi_;
    dg::ClonePtr<dg::aGeometry2d> m_g;
};
 
//used in constructor of Fieldaligned
template<class real_type>
void integrate_all_fieldlines2d( const dg::geo::CylindricalVectorLvl1& vec,
    const dg::aRealGeometry2d<real_type>& grid_field,
    const dg::aRealTopology2d<real_type>& grid_evaluate,
    std::array<thrust::host_vector<real_type>,3>& yp,
    const thrust::host_vector<double>& vol0,
    thrust::host_vector<real_type>& yp2b,
    thrust::host_vector<bool>& in_boxp,
    real_type deltaPhi, real_type eps)
{
    //grid_field contains the global geometry for the field and the boundaries
    //grid_evaluate contains the points to actually integrate
    std::array<thrust::host_vector<real_type>,3> y{
        dg::evaluate( dg::cooX2d, grid_evaluate),
        dg::evaluate( dg::cooY2d, grid_evaluate),
        vol0
    };
    yp.fill(dg::evaluate( dg::zero, grid_evaluate));
    //construct field on high polynomial grid, then integrate it
    dg::geo::detail::DSField field;
    if( !dynamic_cast<const dg::CartesianGrid2d*>( &grid_field))
        field = dg::geo::detail::DSField( vec, grid_field);
 
    //field in case of cartesian grid
    dg::geo::detail::DSFieldCylindrical4 cyl_field(vec);
    const unsigned size = grid_evaluate.size();
    dg::Adaptive<dg::ERKStep<std::array<real_type,3>>> adapt(
            "Dormand-Prince-7-4-5", std::array<real_type,3>{0,0,0});
    dg::AdaptiveTimeloop< std::array<real_type,3>> odeint;
    if( dynamic_cast<const dg::CartesianGrid2d*>( &grid_field))
        odeint = dg::AdaptiveTimeloop<std::array<real_type,3>>( adapt,
                cyl_field, dg::pid_control, dg::fast_l2norm, eps, 1e-10);
    else
        odeint = dg::AdaptiveTimeloop<std::array<real_type,3>>( adapt,
                field, dg::pid_control, dg::fast_l2norm, eps, 1e-10);
 
    for( unsigned i=0; i<size; i++)
    {
        std::array<real_type,3> coords{y[0][i],y[1][i],y[2][i]}, coordsP;
        //x,y,s
        real_type phi1 = deltaPhi;
        odeint.set_dt( deltaPhi/2.);
        odeint.integrate( 0, coords, phi1, coordsP);
        yp[0][i] = coordsP[0], yp[1][i] = coordsP[1], yp[2][i] = coordsP[2];
    }
    yp2b.assign( grid_evaluate.size(), deltaPhi); //allocate memory for output
    in_boxp.resize( yp2b.size());
    //Now integrate again but this time find the boundary distance
    for( unsigned i=0; i<size; i++)
    {
        std::array<real_type,3> coords{y[0][i],y[1][i],y[2][i]}, coordsP;
        in_boxp[i] = grid_field.contains( std::array{yp[0][i], yp[1][i]}) ? true : false;
        if( false == in_boxp[i])
        {
            //x,y,s
            real_type phi1 = deltaPhi;
            odeint.integrate_in_domain( 0., coords, phi1, coordsP, 0., (const
                        dg::aRealTopology2d<real_type>&)grid_field, eps);
            yp2b[i] = phi1;
        }
    }
}
 
// Add non-overlapping matrices from row blocks
template<class real_type>
void add_rowwise_from_sub( dg::IHMatrix_t<real_type>& mat, const dg::IHMatrix_t<real_type>& sub)
{
    // all matrices need same number of cols
    if( sub.num_rows() != mat.num_rows())
        throw dg::Error(dg::Message(_ping_)<<" Submatrix has "<<sub.num_rows()<<" but should have "<<mat.num_rows()<<" columns!");
    if( sub.num_cols() != mat.num_cols())
        throw dg::Error(dg::Message(_ping_)<<" Submatrix has "<<sub.num_cols()<<" but should have "<<mat.num_cols()<<" columns!");
    dg::blas1::axpby( 1, sub.row_offsets(), 1, mat.row_offsets());
    mat.column_indices().insert( mat.column_indices().end(),
        sub.column_indices().begin(), sub.column_indices().end());
    mat.values().insert( mat.values().end(),
        sub.values().begin(), sub.values().end());
}
// Concatenate matrices from row blocks
template<class real_type>
void add_from_sub( dg::IHMatrix_t<real_type>& mat, const dg::IHMatrix_t<real_type>& sub, std::string method)
{
    if( mat.values().empty())
        mat = sub;
    else if( method != "linear")
        add_rowwise_from_sub( mat, sub);
    else
        mat += sub;
    return;
}
 
}//namespace detail
 
 
template<class ProductGeometry, class IMatrix, class container >
struct Fieldaligned
{
 
    Fieldaligned(){}
    template <class Limiter>
    Fieldaligned(const dg::geo::TokamakMagneticField& vec,
        const ProductGeometry& grid,
        dg::bc bcx = dg::NEU,
        dg::bc bcy = dg::NEU,
        Limiter limit = FullLimiter(),
        double eps = 1e-5,
        unsigned mx=12, unsigned my=12,
        double deltaPhi = -1,
        std::string interpolation_method = "linear-nearest",
        bool benchmark=true
        ):
            Fieldaligned( dg::geo::createBHat(vec),
                grid, bcx, bcy, limit, eps, mx, my, deltaPhi, interpolation_method, benchmark)
    {
    }
 
    template <class Limiter>
    Fieldaligned(const dg::geo::CylindricalVectorLvl1& vec,
        const ProductGeometry& grid,
        dg::bc bcx = dg::NEU,
        dg::bc bcy = dg::NEU,
        Limiter limit = FullLimiter(),
        double eps = 1e-5,
        unsigned mx=12, unsigned my=12,
        double deltaPhi = -1,
        std::string interpolation_method = "linear-nearest",
        bool benchmark=true
        );
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = Fieldaligned( std::forward<Params>( ps)...);
    }
 
    dg::bc bcx()const{
        return m_bcx;
    }
    dg::bc bcy()const{
        return m_bcy;
    }
 
 
    void set_boundaries( dg::bc bcz, double left, double right)
    {
        m_bcz = bcz;
        const dg::Grid1d g2d( 0, 1, 1, m_perp_size);
        m_left  = dg::evaluate( dg::CONSTANT(left), g2d);
        m_right = dg::evaluate( dg::CONSTANT(right),g2d);
    }
 
    void set_boundaries( dg::bc bcz, const container& left, const container& right)
    {
        m_bcz = bcz;
        m_left = left;
        m_right = right;
    }
 
    void set_boundaries( dg::bc bcz, const container& global, double scal_left, double scal_right)
    {
        dg::split( global, m_split_in, *m_g);
        dg::blas1::axpby( scal_left,  m_split_in[0],      0, m_left);
        dg::blas1::axpby( scal_right, m_split_in[m_Nz-1], 0, m_right);
        m_bcz = bcz;
    }
 
    void operator()(enum whichMatrix which, const container& in, container& out) const;
 
    double deltaPhi() const{return m_deltaPhi;}
    const container& hbm()const {
        return m_hbm;
    }
    const container& hbp()const {
        return m_hbp;
    }
    const container& sqrtG()const {
        return m_G;
    }
    const container& sqrtGm()const {
        return m_Gm;
    }
    const container& sqrtGp()const {
        return m_Gp;
    }
    const container& bphi()const {
        return m_bphi;
    }
    const container& bphiM()const {
        return m_bphiM;
    }
    const container& bphiP()const {
        return m_bphiP;
    }
    const container& bbm()const {
        return m_bbm;
    }
    const container& bbo()const {
        return m_bbo;
    }
    const container& bbp()const {
        return m_bbp;
    }
    const ProductGeometry& grid()const{return *m_g;}
 
    container interpolate_from_coarse_grid( const ProductGeometry& grid_coarse, const container& coarse) const;
    void integrate_between_coarse_grid( const ProductGeometry& grid_coarse, const container& coarse, container& out ) const;
 
 
    template< class BinaryOp, class UnaryOp>
    container evaluate( BinaryOp binary, UnaryOp unary,
            unsigned p0, unsigned rounds) const;
 
    std::string method() const{return m_interpolation_method;}
 
    private:
    void ePlus( enum whichMatrix which, const container& in, container& out) const;
    void eMinus(enum whichMatrix which, const container& in, container& out) const;
    void zero( enum whichMatrix which, const container& in, container& out) const;
    IMatrix m_plus, m_zero, m_minus; //2d interpolation matrices
    mutable IMatrix m_plusT, m_minusT; // only allocated if necessary
    // backproject
    typename detail::Container2EllSparseBlockMat<container>::type m_bx, m_by, m_inv_bx, m_inv_by;
    bool apply_backproject = false, apply_inv_backproject = false;
    mutable container m_temp0, m_temp1; // 3d size
    //
    container m_hbm, m_hbp;         //3d size
    container m_G, m_Gm, m_Gp; // 3d size
    container m_bphi, m_bphiM, m_bphiP; // 3d size
    container m_bbm, m_bbp, m_bbo;  //3d size masks
 
    container m_left, m_right;      //perp_size
    container m_limiter;            //perp_size
    mutable container m_ghostM, m_ghostP;   //perp_size
    unsigned m_Nz, m_perp_size;
    dg::bc m_bcx, m_bcy, m_bcz;
    mutable std::vector<dg::View<const container>> m_split_in;
    mutable std::vector<dg::View< container>> m_split_out;
    dg::ClonePtr<ProductGeometry> m_g;
    dg::InverseKroneckerTriDiagonal2d<container> m_inv_linear;
    double m_deltaPhi;
    std::string m_interpolation_method;
 
    mutable bool m_have_adjoint = false;
    void updateAdjoint( ) const // only changes mutable m_plusT, m_minusT, m_have_adjoint
    {
        m_plusT = m_plus.transpose();
        m_minusT = m_minus.transpose();
        m_have_adjoint = true;
    }
};
 
 
template<class Geometry, class IMatrix, class container>
template <class Limiter>
Fieldaligned<Geometry, IMatrix, container>::Fieldaligned(
    const dg::geo::CylindricalVectorLvl1& vec,
    const Geometry& grid,
    dg::bc bcx, dg::bc bcy, Limiter limit, double eps,
    unsigned mx, unsigned my, double deltaPhi, std::string interpolation_method, bool benchmark) :
        m_g(grid),
        m_interpolation_method(interpolation_method)
{
 
    std::string inter_m, project_m, fine_m;
    detail::parse_method( interpolation_method, inter_m, project_m, fine_m);
    // For project method "const" we round up to the nearest multiple of n
    if( project_m != "dg" && fine_m == "dg")
    {
        unsigned rx = mx % grid.nx(), ry = my % grid.ny();
        if( 0 != rx || 0 != ry)
        {
            std::cerr << "#Warning: for projection method \"const\" mx and my "
            <<mx<<" "<<my<<" must be multiples of nx and ny "
            <<grid.nx()<<" "<<grid.ny()<<" ! Rounding up for you ...\n";
            mx = mx + grid.nx() - rx;
            my = my + grid.ny() - ry;
        }
    }
    if( benchmark)
        std::cout << "# Interpolation method: \""<<inter_m
            << "\" projection method: \""<<project_m
            <<"\" fine grid \""<<fine_m<<"\"\n";
    if( (grid.bcx() == PER && bcx != PER) || (grid.bcx() != PER && bcx == PER) )
        throw( dg::Error(dg::Message(_ping_)<<"Fieldaligned: Got conflicting periodicity in x. The grid says "<<bc2str(grid.bcx())<<" while the parameter says "<<bc2str(bcx)));
    if( (grid.bcy() == PER && bcy != PER) || (grid.bcy() != PER && bcy == PER) )
        throw( dg::Error(dg::Message(_ping_)<<"Fieldaligned: Got conflicting boundary conditions in y. The grid says "<<bc2str(grid.bcy())<<" while the parameter says "<<bc2str(bcy)));
    m_Nz=grid.Nz(), m_bcx = bcx, m_bcy = bcy, m_bcz=grid.bcz();
    if( deltaPhi <=0) deltaPhi = grid.hz();
    //  grid_trafo -> grid_equi -> grid_fine -> grid_equi -> grid_trafo
    dg::Timer t;
    if( benchmark) t.tic();
    dg::ClonePtr<dg::aGeometry2d> grid_transform( grid.perp_grid()) ;
    std::array<thrust::host_vector<double>,3> yp_trafo, ym_trafo;
    thrust::host_vector<bool> in_boxp, in_boxm;
    thrust::host_vector<double> hbp, hbm;
    thrust::host_vector<double> vol = dg::tensor::volume(grid.metric()), vol2d0;
    auto vol2d = dg::split( vol, grid);
    dg::assign( vol2d[0], vol2d0);
    dg::ClonePtr<dg::aGeometry2d> grid_magnetic = grid_transform;//INTEGRATE HIGH ORDER GRID
    // grid_magnetic is only used for integrating in curvilinear coords
    grid_magnetic->set( grid_transform->n() < 3 ? 4 : 7, grid_magnetic->Nx(), grid_magnetic->Ny());
    detail::integrate_all_fieldlines2d( vec, *grid_magnetic, *grid_transform,
            yp_trafo, vol2d0, hbp, in_boxp, deltaPhi, eps);
    detail::integrate_all_fieldlines2d( vec, *grid_magnetic, *grid_transform,
            ym_trafo, vol2d0, hbm, in_boxm, -deltaPhi, eps);
    if( benchmark)
    {
        t.toc();
        std::cout << "# DS: Fieldline integration took: "<<t.diff()<<"\n";
        t.tic();
    }
    // Assemble minus, zero and plus through sub-grids
    // The idea for the sub-grids is that the fine interpolation matrix takes
    // up a large chunk of memory which can be avoided by sub-dividing the
    // fine grid along its rows
    dg::IHMatrix plus, zero, minus;
    dg::IHMatrix interpolate, zero_interpolate, projection;
    for( unsigned sub = 0; sub < grid_transform->Ny(); sub++)
    {
        dg::RealGrid2d<double> grid_fine_sub(
            grid_transform->x0(),
            grid_transform->x1(),
            grid_transform->y0() + sub*grid_transform->hy(),
            // Fix: construction bug The local right boundary should be the same as the global right boundary
            sub == grid_transform->Ny()-1 ? grid_transform->y1() :
                grid_transform->y0() + (sub+1)*grid_transform->hy(),
            grid_transform->n(), grid_transform->Nx(), 1,
            grid_transform->bcx(), grid_transform->bcy());
        // We do not need metric of grid_equidist or of grid_fine
        dg::RealGrid2d<double> grid_equidist( *grid_transform) ;
        grid_equidist.set( 1, grid_transform->shape(0), grid_transform->shape(1));
        if( fine_m == "equi")
            grid_fine_sub.set( 1, grid_fine_sub.shape(0), grid_fine_sub.shape(1));
        grid_fine_sub.multiplyCellNumbers((double)mx, (double)my);
        std::array<thrust::host_vector<double>,3> yp, ym;
        dg::HVec Xf = dg::evaluate(  dg::cooX2d, grid_fine_sub);
        dg::HVec Yf = dg::evaluate(  dg::cooY2d, grid_fine_sub);
        // interpolate matrix is same in all sub the rows just shift ...
        unsigned shift = grid_transform->shape(0) * grid_transform->n();
        if( sub <= 1 || sub >= grid_transform->Ny() - 2 || grid_transform->n() < 3) // in latter case boundary conditions could destroy invariance
        {
            interpolate = dg::create::interpolation( Xf, Yf,
                *grid_transform, dg::NEU, dg::NEU, grid_transform->n() < 3 ? "cubic" : "dg");
            zero_interpolate = dg::create::interpolation( Xf, Yf,
                inter_m == "dg" ? dg::RealGrid2d<double>(*grid_transform) :
                grid_equidist, bcx, bcy, inter_m);
        }
        else
        {
            dg::blas1::plus( interpolate.column_indices(), shift);
            dg::blas1::plus( zero_interpolate.column_indices(), shift);
        }
        yp.fill(dg::evaluate( dg::zero, grid_fine_sub));
        ym = yp;
        for( int i=0; i<2; i++) //only R and Z get interpolated
        {
            dg::blas2::symv( interpolate, yp_trafo[i], yp[i]);
            dg::blas2::symv( interpolate, ym_trafo[i], ym[i]);
        }
        if( sub <= 1 || sub >= grid_transform->Ny() -2 )
        {
            if( project_m == "dg")
            {
                // Note that it is possible to project onto a bigger grid, the corresponding rows are just empty
                projection = dg::create::projection( *grid_transform, grid_fine_sub);
            }
            else
            {
                projection = dg::create::projection( grid_equidist,
                    grid_fine_sub, project_m);
            }
        }
        else
        {
            // how to add to rows in csr formatted matrix:
            projection.row_offsets().insert( projection.row_offsets().begin(),
                shift, 0);
            projection.row_offsets().erase( projection.row_offsets().end() -
                shift, projection.row_offsets().end());
        }
        std::array<dg::HVec*,3> xcomp{ &yp[0], &Xf, &ym[0]};
        std::array<dg::HVec*,3> ycomp{ &yp[1], &Yf, &ym[1]};
        std::array<dg::IHMatrix*,3> result{ &plus, &zero, &minus};
        dg::IHMatrix subresult;
 
        for( unsigned u=0; u<3; u++)
        {
            if( u == 1)
                subresult = projection*zero_interpolate;
            else
                subresult = projection*dg::create::interpolation( *xcomp[u], *ycomp[u],
                    inter_m == "dg" ? dg::RealGrid2d<double>(*grid_transform) :
                    grid_equidist, bcx, bcy, inter_m);
            detail::add_from_sub( *result[u], subresult, project_m);
        }
    }
    // 26.9.25: It turns out that directly sandwiching the interpolation
    // matrices minus,zero,plus with B^{-1} I B makes the average number of
    // non-zeros per row explode (for n=3 from 9 to up to 50 - 75!). Therefore
    // it is faster to separately store B^{-1} and B and apply them before and
    // after I. Separate benchmarks of B = Bx x By reveals that 2 of our
    // sparseblockmat symv implementation (Bx and By separately) is faster than
    // any attempt at merging the two kernels into one and certainly faster
    // than the corresponding dg::IDMatrix symv kernel. In fact, the 2 symv
    // kernels are almost as fast as just one axpby. So, here we go:
    if( project_m != "dg")
    {
        apply_inv_backproject = true; // from equidist to dg
        m_inv_bx = dg::create::fast_transform( 0,
            dg::invert(dg::create::detail::square_backproject( grid.grid(0))),
            grid);
        m_inv_by = dg::create::fast_transform( 1,
            dg::invert(dg::create::detail::square_backproject( grid.grid(1))),
            grid);
    }
    if( inter_m != "dg")
    {
        apply_backproject = true; // from dg to equidist
        m_bx = dg::create::fast_transform( 0,
            dg::create::detail::square_backproject( grid.grid(0)),
            grid);
        m_by = dg::create::fast_transform( 1,
            dg::create::detail::square_backproject( grid.grid(1)),
            grid);
    }
    dg::blas2::transfer( minus, m_minus);
    dg::blas2::transfer( zero, m_zero);
    dg::blas2::transfer( plus, m_plus);
    m_have_adjoint = false;
 
    if( benchmark)
    {
        t.toc();
        std::cout << "# DS: Assembly of matrices  took: "<<t.diff()<<"\n";
        std::streamsize ss = std::cout.precision();
        std::cout << std::setprecision(1) << std::fixed;
        std::cout << "# DS: Average nnz per row (plus | zero | minus): "
                  //<<(double)backproject.num_nnz()/(double)backproject.num_rows()<<" | "
                  <<(double)plus.num_nnz()/(double)plus.num_rows()<<" | "
                  <<(double)zero.num_nnz()/(double)zero.num_rows()<<" | "
                  <<(double)minus.num_nnz()/(double)minus.num_rows()<<"\n";
        std::cout << std::defaultfloat << std::setprecision(ss);
    }
    dg::HVec hbphi( yp_trafo[2]), hbphiP(hbphi), hbphiM(hbphi);
    hbphi = dg::pullback( vec.z(), *grid_transform);
    //this is a pullback bphi( R(zeta, eta), Z(zeta, eta)):
    if( dynamic_cast<const dg::CartesianGrid2d*>( grid_transform.get()))
    {
        for( unsigned i=0; i<hbphiP.size(); i++)
        {
            hbphiP[i] = vec.z()(yp_trafo[0][i], yp_trafo[1][i]);
            hbphiM[i] = vec.z()(ym_trafo[0][i], ym_trafo[1][i]);
        }
    }
    else
    {
        dg::HVec Ihbphi = dg::pullback( vec.z(), *grid_magnetic);
        dg::HVec Lhbphi = dg::forward_transform( Ihbphi, *grid_magnetic);
        for( unsigned i=0; i<yp_trafo[0].size(); i++)
        {
            hbphiP[i] = dg::interpolate( dg::lspace, Lhbphi, yp_trafo[0][i],
                    yp_trafo[1][i], *grid_magnetic);
            hbphiM[i] = dg::interpolate( dg::lspace, Lhbphi, ym_trafo[0][i],
                    ym_trafo[1][i], *grid_magnetic);
        }
    }
    dg::assign3dfrom2d( hbphi,  m_bphi,  grid);
    dg::assign3dfrom2d( hbphiM, m_bphiM, grid);
    dg::assign3dfrom2d( hbphiP, m_bphiP, grid);
 
    dg::assign3dfrom2d( yp_trafo[2], m_Gp, grid);
    dg::assign3dfrom2d( ym_trafo[2], m_Gm, grid);
    // The weights don't matter since they fall out in Div and Lap anyway
    // But they are good for testing
    m_G = vol;
    container weights = dg::create::weights( grid);
    dg::blas1::pointwiseDot( m_G, weights, m_G);
    dg::blas1::pointwiseDot( m_Gp, weights, m_Gp);
    dg::blas1::pointwiseDot( m_Gm, weights, m_Gm);
 
    dg::assign( dg::evaluate( dg::zero, grid), m_hbm);
    m_split_in     = dg::split( (const container&)m_hbm, grid);
    m_split_out  = dg::split( m_hbm, grid);
    dg::assign3dfrom2d( hbp, m_hbp, grid);
    dg::assign3dfrom2d( hbm, m_hbm, grid);
    dg::blas1::scal( m_hbm, -1.);
 
    m_temp0 = m_temp1 = m_hbm; // 3d size
 
    thrust::host_vector<double> bbm( in_boxp.size(),0.), bbo(bbm), bbp(bbm);
    for( unsigned i=0; i<in_boxp.size(); i++)
    {
        if( !in_boxp[i] && !in_boxm[i])
            bbo[i] = 1.;
        else if( !in_boxp[i] && in_boxm[i])
            bbp[i] = 1.;
        else if( in_boxp[i] && !in_boxm[i])
            bbm[i] = 1.;
        // else all are 0
    }
    dg::assign3dfrom2d( bbm, m_bbm, grid);
    dg::assign3dfrom2d( bbo, m_bbo, grid);
    dg::assign3dfrom2d( bbp, m_bbp, grid);
 
    m_deltaPhi = deltaPhi; // store for evaluate
 
    m_perp_size = grid_transform->size();
    dg::assign( dg::pullback(limit, *grid_transform), m_limiter);
    dg::assign( dg::evaluate(dg::zero, *grid_transform), m_left);
    m_ghostM = m_ghostP = m_right = m_left;
}
 
 
template<class G, class I, class container>
container Fieldaligned<G, I,container>::interpolate_from_coarse_grid(
        const G& grid, const container& in) const
{
    //I think we need grid as input to split input vector and we need to interpret
    //the grid nodes as node centered not cell-centered!
    //idea: apply I+/I- cphi - 1 times in each direction and then apply interpolation formula
    assert( m_g->Nz() % grid.Nz() == 0);
    unsigned Nz_coarse = grid.Nz(), Nz = m_g->Nz();
    unsigned cphi = Nz / Nz_coarse;
 
    container out = dg::evaluate( dg::zero, *m_g);
    container helper = dg::evaluate( dg::zero, *m_g);
    dg::split( helper, m_split_out, *m_g);
    std::vector<dg::View< container>> out_split = dg::split( out, *m_g);
    std::vector<dg::View< const container>> in_split = dg::split( in, grid);
    if( apply_backproject)
    {
        dg::blas2::symv( m_bx, in, m_temp0);
        dg::blas2::symv( m_by, m_temp0, m_temp1);
        in_split = dg::split( m_temp1, grid);
    }
    for ( int i=0; i<(int)Nz_coarse; i++)
    {
        //1. copy input vector to appropriate place in output
        dg::blas1::copy( in_split[i], out_split[i*cphi]);
        dg::blas1::copy( in_split[i], m_split_out[i*cphi]);
    }
    //Step 1 needs to finish so that m_split_out contains values everywhere
    //2. Now apply plus and minus T to fill in the rest
    for ( int i=0; i<(int)Nz_coarse; i++)
    {
        for( int j=1; j<(int)cphi; j++)
        {
            dg::blas2::symv( m_minus, out_split[i*cphi+j-1], out_split[i*cphi+j]);
            dg::blas2::symv( m_plus, m_split_out[(i*cphi+cphi+1-j)%Nz], m_split_out[i*cphi+cphi-j]);
        }
    }
    //3. Now add up with appropriate weights
    for( int i=0; i<(int)Nz_coarse; i++)
        for( int j=1; j<(int)cphi; j++)
        {
            double alpha = (double)(cphi-j)/(double)cphi;
            double beta = (double)j/(double)cphi;
            dg::blas1::axpby( alpha, out_split[i*cphi+j], beta, m_split_out[i*cphi+j], out_split[i*cphi+j]);
        }
    if( apply_inv_backproject)
    {
        dg::blas2::symv( m_inv_bx, out, m_temp0);
        dg::blas2::symv( m_inv_by, m_temp0, out);
    }
    return out;
}
template<class G, class I, class container>
void Fieldaligned<G, I,container>::integrate_between_coarse_grid( const G& grid, const container& in, container& out) const
{
    // TODO: This function is missing a test
    assert( m_g->Nz() % grid.Nz() == 0);
    unsigned Nz_coarse = grid.Nz(), Nz = m_g->Nz();
    unsigned cphi = Nz / Nz_coarse;
 
    out = in;
    container helperP( in), helperM(in), tempP(in), tempM(in);
    typename detail::Container2EllSparseBlockMat<container>::type bx, by, inv_bx, inv_by;
    if( apply_backproject)
    {
        bx = dg::create::fast_transform( 0,
            dg::create::detail::square_backproject( grid.grid(0)),
            grid);
        by = dg::create::fast_transform( 1,
            dg::create::detail::square_backproject( grid.grid(1)),
            grid);
    }
    if( apply_inv_backproject)
    {
        inv_bx = dg::create::fast_transform( 0,
            dg::invert(dg::create::detail::square_backproject( grid.grid(0))),
            grid);
        inv_by = dg::create::fast_transform( 1,
            dg::invert(dg::create::detail::square_backproject( grid.grid(1))),
            grid);
    }
 
    //1. Apply plus and minus T and sum up
    for( int j=1; j<(int)cphi; j++)
    {
        if( apply_backproject)
        {
            dg::blas2::symv( bx, helperP, tempP);
            dg::blas2::symv( by, tempP, helperP);
            dg::blas2::symv( bx, helperM, tempM);
            dg::blas2::symv( by, tempM, helperM);
        }
        dg::blas2::symv( m_minus, helperP, tempP);
        if( apply_inv_backproject)
        {
            dg::blas2::symv( inv_bx, tempP, helperP);
            dg::blas2::symv( inv_by, helperP, tempP);
        }
        dg::blas1::axpby( (double)(cphi-j)/(double)cphi, tempP, 1., out  );
        helperP.swap(tempP);
        dg::blas2::symv( m_plus, helperM, tempM);
        if( apply_inv_backproject)
        {
            dg::blas2::symv( inv_bx, tempM, helperM);
            dg::blas2::symv( inv_by, helperM, tempM);
        }
        dg::blas1::axpby( (double)(cphi-j)/(double)cphi, tempM, 1., out  );
        helperM.swap(tempM);
    }
    dg::blas1::scal( out, 1./(double)cphi);
}
 
template<class G, class I, class container>
void Fieldaligned<G, I, container >::operator()(enum whichMatrix which, const container& f, container& fe) const
{
    if(     which == einsPlus  || which == einsMinusT ) ePlus(  which, f, fe);
    else if(which == einsMinus || which == einsPlusT  ) eMinus( which, f, fe);
    else if(which == zeroMinus || which == zeroPlus ||
            which == zeroMinusT|| which == zeroPlusT ||
            which == zeroForw  ) zero(   which, f, fe);
}
 
template< class G, class I, class container>
void Fieldaligned<G, I, container>::zero( enum whichMatrix which,
        const container& f, container& f0) const
{
    dg::split( f, m_split_in, *m_g);
    dg::split( f0, m_split_out, *m_g);
    if( apply_backproject)
    {
        dg::blas2::symv( m_bx, f, m_temp0);
        dg::blas2::symv( m_by, m_temp0, m_temp1);
        dg::split( (const container&)m_temp1, m_split_in, *m_g);
    }
    //1. compute 2d interpolation in every plane and store in m_split_out
    for( unsigned i0=0; i0<m_Nz; i0++)
    {
        if(which == zeroPlus)
            dg::blas2::symv( m_plus,   m_split_in[i0], m_split_out[i0]);
        else if(which == zeroMinus)
            dg::blas2::symv( m_minus,  m_split_in[i0], m_split_out[i0]);
        else if(which == zeroPlusT)
        {
            if( ! m_have_adjoint) updateAdjoint( );
            dg::blas2::symv( m_plusT,  m_split_in[i0], m_split_out[i0]);
        }
        else if(which == zeroMinusT)
        {
            if( ! m_have_adjoint) updateAdjoint( );
            dg::blas2::symv( m_minusT, m_split_in[i0], m_split_out[i0]);
        }
        else if( which == zeroForw)
        {
            if ( m_interpolation_method != "dg" )
            {
                dg::blas2::symv( m_zero, m_split_in[i0], m_split_out[i0]);
            }
            else
                dg::blas1::copy( m_split_in[i0], m_split_out[i0]);
        }
    }
    if( apply_inv_backproject)
    {
        dg::blas2::symv( m_inv_bx, f0, m_temp0);
        dg::blas2::symv( m_inv_by, m_temp0, f0);
    }
}
template< class G, class I, class container>
void Fieldaligned<G, I, container>::ePlus( enum whichMatrix which,
        const container& f, container& fpe) const
{
    dg::split( f, m_split_in, *m_g);
    dg::split( fpe, m_split_out, *m_g);
    if( apply_backproject)
    {
        dg::blas2::symv( m_bx, f, m_temp0);
        dg::blas2::symv( m_by, m_temp0, m_temp1);
        dg::split( (const container&)m_temp1, m_split_in, *m_g);
    }
    //1. compute 2d interpolation in every plane and store in m_split_out
    for( unsigned i0=0; i0<m_Nz; i0++)
    {
        unsigned ip = (i0==m_Nz-1) ? 0:i0+1;
        if(which == einsPlus)
            dg::blas2::symv( m_plus,   m_split_in[ip], m_split_out[i0]);
        else if(which == einsMinusT)
        {
            if( ! m_have_adjoint) updateAdjoint( );
            dg::blas2::symv( m_minusT, m_split_in[ip], m_split_out[i0]);
        }
    }
    if( apply_inv_backproject)
    {
        dg::blas2::symv( m_inv_bx, fpe, m_temp0);
        dg::blas2::symv( m_inv_by, m_temp0, fpe);
    }
    //2. apply right boundary conditions in last plane
    unsigned i0=m_Nz-1;
    if( m_bcz != dg::PER)
    {
        dg::split( f, m_split_in, *m_g);
        if( m_bcz == dg::DIR || m_bcz == dg::NEU_DIR)
            dg::blas1::axpby( 2, m_right, -1., m_split_in[i0], m_ghostP);
        if( m_bcz == dg::NEU || m_bcz == dg::DIR_NEU)
            dg::blas1::axpby( m_deltaPhi, m_right, 1., m_split_in[i0], m_ghostP);
        //interlay ghostcells with periodic cells: L*g + (1-L)*fpe
        dg::blas1::axpby( 1., m_ghostP, -1., m_split_out[i0], m_ghostP);
        dg::blas1::pointwiseDot( 1., m_limiter, m_ghostP, 1., m_split_out[i0]);
    }
}
 
template< class G, class I, class container>
void Fieldaligned<G, I, container>::eMinus( enum whichMatrix which,
        const container& f, container& fme) const
{
    dg::split( f, m_split_in, *m_g);
    dg::split( fme, m_split_out, *m_g);
    if( apply_backproject)
    {
        dg::blas2::symv( m_bx, f, m_temp0);
        dg::blas2::symv( m_by, m_temp0, m_temp1);
        dg::split( (const container&)m_temp1, m_split_in, *m_g);
    }
    //1. compute 2d interpolation in every plane and store in m_split_out
    for( unsigned i0=0; i0<m_Nz; i0++)
    {
        unsigned im = (i0==0) ? m_Nz-1:i0-1;
        if(which == einsPlusT)
        {
            if( ! m_have_adjoint) updateAdjoint( );
            dg::blas2::symv( m_plusT, m_split_in[im], m_split_out[i0]);
        }
        else if (which == einsMinus)
            dg::blas2::symv( m_minus, m_split_in[im], m_split_out[i0]);
    }
    if( apply_inv_backproject)
    {
        dg::blas2::symv( m_inv_bx, fme, m_temp0);
        dg::blas2::symv( m_inv_by, m_temp0, fme);
    }
    //2. apply left boundary conditions in first plane
    unsigned i0=0;
    if( m_bcz != dg::PER)
    {
        if( m_bcz == dg::DIR || m_bcz == dg::DIR_NEU)
            dg::blas1::axpby( 2., m_left,  -1., m_split_in[i0], m_ghostM);
        if( m_bcz == dg::NEU || m_bcz == dg::NEU_DIR)
            dg::blas1::axpby( -m_deltaPhi, m_left, 1., m_split_in[i0], m_ghostM);
        //interlay ghostcells with periodic cells: L*g + (1-L)*fme
        dg::blas1::axpby( 1., m_ghostM, -1., m_split_out[i0], m_ghostM);
        dg::blas1::pointwiseDot( 1., m_limiter, m_ghostM, 1., m_split_out[i0]);
    }
}
 
template<class G, class I, class container>
template< class BinaryOp, class UnaryOp>
container Fieldaligned<G, I,container>::evaluate( BinaryOp binary,
        UnaryOp unary, unsigned p0, unsigned rounds) const
{
    //idea: simply apply I+/I- enough times on the init2d vector to get the result in each plane
    //unary function is always such that the p0 plane is at x=0
    assert( p0 < m_g->Nz());
    const dg::ClonePtr<aGeometry2d> g2d = m_g->perp_grid();
    container init2d = dg::pullback( binary, *g2d);
    container zero2d = dg::evaluate( dg::zero, *g2d);
    typename detail::Container2EllSparseBlockMat<container>::type bx, by, inv_bx, inv_by;
    if( apply_backproject)
    {
        bx = dg::create::fast_transform( 0,
            dg::create::detail::square_backproject( g2d->grid(0)),
            *g2d);
        by = dg::create::fast_transform( 1,
            dg::create::detail::square_backproject( g2d->grid(1)),
            *g2d);
    }
    if( apply_inv_backproject)
    {
        inv_bx = dg::create::fast_transform( 0,
            dg::invert(dg::create::detail::square_backproject( g2d->grid(0))),
            *g2d);
        inv_by = dg::create::fast_transform( 1,
            dg::invert(dg::create::detail::square_backproject( g2d->grid(1))),
            *g2d);
    }
 
    container temp(init2d), tempP(init2d), tempM(init2d);
    container vec3d = dg::evaluate( dg::zero, *m_g);
    std::vector<container>  plus2d(m_Nz, zero2d), minus2d(plus2d), result(plus2d);
    unsigned turns = rounds;
    if( turns ==0) turns++;
    //first apply Interpolation many times, scale and store results
    for( unsigned r=0; r<turns; r++)
        for( unsigned i0=0; i0<m_Nz; i0++)
        {
            dg::blas1::copy( init2d, tempP);
            dg::blas1::copy( init2d, tempM);
            unsigned rep = r*m_Nz + i0;
            for(unsigned k=0; k<rep; k++)
            {
                if( apply_backproject)
                {
                    dg::blas2::symv( bx, tempP, temp);
                    dg::blas2::symv( by, temp, tempP);
                    dg::blas2::symv( bx, tempM, temp);
                    dg::blas2::symv( by, temp, tempM);
                }
                dg::blas2::symv( m_minus, tempP, temp);
                if( apply_inv_backproject)
                {
                    dg::blas2::symv( inv_bx, temp, tempP);
                    dg::blas2::symv( inv_by, tempP, temp);
                }
                temp.swap( tempP);
                dg::blas2::symv( m_plus, tempM, temp);
                if( apply_inv_backproject)
                {
                    dg::blas2::symv( inv_bx, temp, tempM);
                    dg::blas2::symv( inv_by, tempM, temp);
                }
                temp.swap( tempM);
            }
            dg::blas1::scal( tempP, unary(  (double)rep*m_deltaPhi ) );
            dg::blas1::scal( tempM, unary( -(double)rep*m_deltaPhi ) );
            dg::blas1::axpby( 1., tempP, 1., plus2d[i0]);
            dg::blas1::axpby( 1., tempM, 1., minus2d[i0]);
        }
    //now we have the plus and the minus filaments
    if( rounds == 0) //there is a limiter
    {
        for( unsigned i0=0; i0<m_Nz; i0++)
        {
            int idx = (int)i0 - (int)p0;
            if(idx>=0)
                result[i0] = plus2d[idx];
            else
                result[i0] = minus2d[abs(idx)];
            thrust::copy( result[i0].begin(), result[i0].end(), vec3d.begin() + i0*m_perp_size);
        }
    }
    else //sum up plus2d and minus2d
    {
        for( unsigned i0=0; i0<m_Nz; i0++)
        {
            unsigned revi0 = (m_Nz - i0)%m_Nz; //reverted index
            dg::blas1::axpby( 1., plus2d[i0], 0., result[i0]);
            dg::blas1::axpby( 1., minus2d[revi0], 1., result[i0]);
        }
        dg::blas1::axpby( -1., init2d, 1., result[0]);
        for(unsigned i0=0; i0<m_Nz; i0++)
        {
            int idx = ((int)i0 -(int)p0 + m_Nz)%m_Nz; //shift index
            thrust::copy( result[idx].begin(), result[idx].end(), vec3d.begin() + i0*m_perp_size);
        }
    }
    return vec3d;
}
 
 
 
template<class BinaryOp, class UnaryOp>
thrust::host_vector<double> fieldaligned_evaluate(
        const aProductGeometry3d& grid,
        const CylindricalVectorLvl0& vec,
        const BinaryOp& binary,
        const UnaryOp& unary,
        unsigned p0,
        unsigned rounds,
        double eps = 1e-5)
{
    unsigned Nz = grid.Nz();
    const dg::ClonePtr<aGeometry2d> g2d = grid.perp_grid();
    // Construct for field-aligned output
    dg::HVec tempP = dg::evaluate( dg::zero, *g2d), tempM( tempP);
    std::vector<dg::HVec>  plus2d(Nz, tempP), minus2d(plus2d), result(plus2d);
    dg::HVec vec3d = dg::evaluate( dg::zero, grid);
    dg::HVec init2d = dg::pullback( binary, *g2d);
    std::array<dg::HVec,3> yy0{
        dg::pullback( dg::cooX2d, *g2d),
        dg::pullback( dg::cooY2d, *g2d),
        dg::evaluate( dg::zero, *g2d)}, yy1(yy0), xx0( yy0), xx1(yy0); //s
    dg::geo::detail::DSFieldCylindrical3 cyl_field(vec);
    double deltaPhi = grid.hz();
    double phiM0 = 0., phiP0 = 0.;
    unsigned turns = rounds;
    if( turns == 0) turns++;
    for( unsigned r=0; r<turns; r++)
        for( unsigned  i0=0; i0<Nz; i0++)
        {
            unsigned rep = r*Nz + i0;
            if( rep == 0)
                tempM = tempP = init2d;
            else
            {
                dg::Adaptive<dg::ERKStep<std::array<double,3>>> adapt(
                        "Dormand-Prince-7-4-5", std::array<double,3>{0,0,0});
                dg::AdaptiveTimeloop<std::array<double,3>> odeint( adapt,
                        cyl_field, dg::pid_control, dg::fast_l2norm, eps,
                        1e-10);
                for( unsigned i=0; i<g2d->size(); i++)
                {
                    // minus direction needs positive integration!
                    double phiM1 = phiM0 + deltaPhi;
                    std::array<double,3>
                        coords0{yy0[0][i],yy0[1][i],yy0[2][i]}, coords1;
                    odeint.integrate_in_domain( phiM0, coords0, phiM1, coords1,
                            deltaPhi, *g2d, eps);
                    yy1[0][i] = coords1[0], yy1[1][i] = coords1[1], yy1[2][i] =
                        coords1[2];
                    tempM[i] = binary( yy1[0][i], yy1[1][i]);
 
                    // plus direction needs negative integration!
                    double phiP1 = phiP0 - deltaPhi;
                    coords0 = std::array<double,3>{xx0[0][i],xx0[1][i],xx0[2][i]};
                    odeint.integrate_in_domain( phiP0, coords0, phiP1, coords1,
                            -deltaPhi, *g2d, eps);
                    xx1[0][i] = coords1[0], xx1[1][i] = coords1[1], xx1[2][i] =
                        coords1[2];
                    tempP[i] = binary( xx1[0][i], xx1[1][i]);
                }
                std::swap( yy0, yy1);
                std::swap( xx0, xx1);
                phiM0 += deltaPhi;
                phiP0 -= deltaPhi;
            }
            dg::blas1::scal( tempM, unary( -(double)rep*deltaPhi ) );
            dg::blas1::scal( tempP, unary(  (double)rep*deltaPhi ) );
            dg::blas1::axpby( 1., tempM, 1., minus2d[i0]);
            dg::blas1::axpby( 1., tempP, 1., plus2d[i0]);
        }
    //now we have the plus and the minus filaments
    if( rounds == 0) //there is a limiter
    {
        for( unsigned i0=0; i0<Nz; i0++)
        {
            int idx = (int)i0 - (int)p0;
            if(idx>=0)
                result[i0] = plus2d[idx];
            else
                result[i0] = minus2d[abs(idx)];
            thrust::copy( result[i0].begin(), result[i0].end(), vec3d.begin() +
                    i0*g2d->size());
        }
    }
    else //sum up plus2d and minus2d
    {
        for( unsigned i0=0; i0<Nz; i0++)
        {
            unsigned revi0 = (Nz - i0)%Nz; //reverted index
            dg::blas1::axpby( 1., plus2d[i0], 0., result[i0]);
            dg::blas1::axpby( 1., minus2d[revi0], 1., result[i0]);
        }
        dg::blas1::axpby( -1., init2d, 1., result[0]);
        for(unsigned i0=0; i0<Nz; i0++)
        {
            int idx = ((int)i0 -(int)p0 + Nz)%Nz; //shift index
            thrust::copy( result[idx].begin(), result[idx].end(), vec3d.begin()
                    + i0*g2d->size());
        }
    }
    return vec3d;
}
 
}//namespace geo
}//namespace dg