FDScheme.hpp

/*
 * FiniteDifferences.hpp
 *
 *  Created on: Sep 17, 2015
 *      Author: i-bird
 */
 
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
 
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Grid/Iterators/grid_dist_id_iterator_sub.hpp"
#include "eq.hpp"
#include "NN/CellList/CellDecomposer.hpp"
#include "Grid/staggered_dist_grid_util.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Vector/Vector_util.hpp"
#include "Grid/staggered_dist_grid.hpp"
 
 
template<typename Sys_eqs>
class FDScheme
{
public:
 
    typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,aggregate<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
 
    typedef Sys_eqs Sys_eqs_typ;
 
    struct constant_b
    {
        typename Sys_eqs::stype scal;
 
        constant_b(typename Sys_eqs::stype scal)
        {
            this->scal = scal;
        }
 
        typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
        {
            return scal;
        }
    };
 
    template<typename grid, unsigned int prp>
    struct grid_b
    {
        grid & gr;
 
        grid_b(grid & gr)
        :gr(gr)
        {}
 
        auto get(grid_dist_key_dx<Sys_eqs::dims> & key) -> decltype(gr.template get<prp>(key))
        {
            return gr.template get<prp>(key);
        }
    };
 
private:
 
    Padding<Sys_eqs::dims> pd;
 
    typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;
 
    typename Sys_eqs::Vector_type b;
 
    const grid_sm<Sys_eqs::dims,void> & gs;
 
    typename Sys_eqs::stype spacing[Sys_eqs::dims];
 
    g_map_type g_map;
 
    size_t row;
 
    size_t row_b;
 
    openfpm::vector<size_t> pnt;
 
    comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar];
 
    size_t s_pnt;
 
    struct key_and_eq
    {
        grid_key_dx<Sys_eqs::dims> key;
 
        size_t eq;
    };
 
    inline key_and_eq from_row_to_key(size_t row)
    {
        key_and_eq ke;
        auto it = g_map.getDomainIterator();
 
        while (it.isNext())
        {
            size_t row_low = g_map.template get<0>(it.get());
 
            if (row >= row_low * Sys_eqs::nvar && row < row_low * Sys_eqs::nvar + Sys_eqs::nvar)
            {
                ke.eq = row - row_low * Sys_eqs::nvar;
                ke.key = g_map.getGKey(it.get());
                return ke;
            }
 
            ++it;
        }
        std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " the row does not map to any position" << "\n";
 
        return ke;
    }
 
    inline const std::vector<size_t> padding( const size_t (& sz)[Sys_eqs::dims], Padding<Sys_eqs::dims> & pd)
    {
        std::vector<size_t> g_sz_pad(Sys_eqs::dims);
 
        for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
            g_sz_pad[i] = sz[i] + pd.getLow(i) + pd.getHigh(i);
 
        return g_sz_pad;
    }
 
    void consistency()
    {
        openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
 
        // A and B must have the same rows
        if (row != row_b)
            std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";
 
        // Indicate all the non zero rows
        openfpm::vector<unsigned char> nz_rows;
        nz_rows.resize(row_b);
 
        for (size_t i = 0 ; i < trpl.size() ; i++)
            nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
 
        // Indicate all the non zero colums
        // This check can be done only on single processor
 
        Vcluster<> & v_cl = create_vcluster();
        if (v_cl.getProcessingUnits() == 1)
        {
            openfpm::vector<unsigned> nz_cols;
            nz_cols.resize(row_b);
 
            for (size_t i = 0 ; i < trpl.size() ; i++)
                nz_cols.get(trpl.get(i).col()) = true;
 
            // all the rows must have a non zero element
            for (size_t i = 0 ; i < nz_rows.size() ; i++)
            {
                if (nz_rows.get(i) == false)
                {
                    key_and_eq ke = from_row_to_key(i);
                    std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i <<  " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
                }
            }
 
            // all the colums must have a non zero element
            for (size_t i = 0 ; i < nz_cols.size() ; i++)
            {
                if (nz_cols.get(i) == false)
                    std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
            }
        }
    }
 
    template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_staggered(Vct & v, Grid_dst & g_dst, grid_key_dx<Grid_dst::dims> & start_d, grid_key_dx<Grid_dst::dims> & stop_d)
    {
        // check that g_dst is staggered
        if (g_dst.is_staggered() == false)
            std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;
 
        // sub-grid iterator over the grid map
        auto g_map_it = g_map.getSubDomainIterator(start_d,stop_d);
 
        // Iterator over the destination grid
        auto g_dst_it = g_dst.getSubDomainIterator(start_d,stop_d);
 
        while (g_map_it.isNext() == true)
        {
            typedef typename to_boost_vmpl<pos...>::type vid;
            typedef boost::mpl::size<vid> v_size;
 
            auto key_src = g_map_it.get();
            size_t lin_id = g_map.template get<0>(key_src);
 
            // destination point
            auto key_dst = g_dst_it.get();
 
            // Transform this id into an id for the Eigen vector
 
            copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());
 
            boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
 
            ++g_map_it;
            ++g_dst_it;
        }
    }
 
    template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_normal(Vct & v, Grid_dst & g_dst)
    {
        // check that g_dst is staggered
        if (g_dst.is_staggered() == true)
            std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be normal " << std::endl;
 
        grid_key_dx<Grid_dst::dims> start;
        grid_key_dx<Grid_dst::dims> stop;
 
        for (size_t i = 0 ; i < Grid_dst::dims ; i++)
        {
            start.set_d(i,pd.getLow(i));
            stop.set_d(i,g_map.size(i) - pd.getHigh(i));
        }
 
        // sub-grid iterator over the grid map
        auto g_map_it = g_map.getSubDomainIterator(start,stop);
 
        // Iterator over the destination grid
        auto g_dst_it = g_dst.getDomainIterator();
 
        while (g_dst_it.isNext() == true)
        {
            typedef typename to_boost_vmpl<pos...>::type vid;
            typedef boost::mpl::size<vid> v_size;
 
            auto key_src = g_map_it.get();
            size_t lin_id = g_map.template get<0>(key_src);
 
            // destination point
            auto key_dst = g_dst_it.get();
 
            // Transform this id into an id for the vector
 
            copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());
 
            boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
 
            ++g_map_it;
            ++g_dst_it;
        }
    }
 
    template<typename bop, typename iterator>
    void impose_dit_b(bop num,
                      long int id ,
                      const iterator & it_d)
    {
 
        auto it = it_d;
        grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
 
        // iterate all the grid points
        while (it.isNext())
        {
            // get the position
            auto key = it.get();
 
            b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
 
            // if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//          T::position(key,gs,s_pos);
#endif
 
            ++row_b;
            ++it;
        }
    }
 
    template<typename T, typename bop, typename iterator> void impose_dit(const T & op,
                                     bop num,
                                     long int id ,
                                     const iterator & it_d)
    {
        openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
 
        auto it = it_d;
        grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
 
        std::unordered_map<long int,typename Sys_eqs::stype> cols;
 
        // iterate all the grid points
        while (it.isNext())
        {
            // get the position
            auto key = it.get();
 
            // Add padding
            for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
                key.getKeyRef().set_d(i,key.getKeyRef().get(i) + pd.getLow(i));
 
            // Calculate the non-zero colums
            T::value(g_map,key,gs,spacing,cols,1.0);
 
            // indicate if the diagonal has been set
            bool is_diag = false;
 
            // create the triplet
            for ( auto it = cols.begin(); it != cols.end(); ++it )
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().col() = it->first;
                trpl.last().value() = it->second;
 
                if (trpl.last().row() == trpl.last().col())
                {is_diag = true;}
 
            }
 
            // If does not have a diagonal entry put it to zero
            if (is_diag == false)
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().value() = 0.0;
            }
 
            b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
 
            cols.clear();
 
            // if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//          T::position(key,gs,s_pos);
#endif
 
            ++row;
            ++row_b;
            ++it;
        }
    }
 
 
 
    void construct_gmap()
    {
        Vcluster<> & v_cl = create_vcluster();
 
        // Calculate the size of the local domain
        size_t sz = g_map.getLocalDomainSize();
 
        // Get the total size of the local grids on each processors
        v_cl.allGather(sz,pnt);
        v_cl.execute();
        s_pnt = 0;
 
        // calculate the starting point for this processor
        for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
            s_pnt += pnt.get(i);
 
        // resize b if needed
        b.resize(Sys_eqs::nvar * g_map.size(),Sys_eqs::nvar * sz);
 
        // Calculate the starting point
 
        // Counter
        size_t cnt = 0;
 
        // Create the re-mapping grid
        auto it = g_map.getDomainIterator();
 
        while (it.isNext())
        {
            auto key = it.get();
 
            g_map.template get<0>(key) = cnt + s_pnt;
 
            ++cnt;
            ++it;
        }
 
        // sync the ghost
        g_map.template ghost_get<0>();
    }
 
    void Initialize(const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain)
    {
        construct_gmap();
 
        // Create a CellDecomposer and calculate the spacing
 
        size_t sz_g[Sys_eqs::dims];
        for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
        {
            if (Sys_eqs::boundary[i] == NON_PERIODIC)
                sz_g[i] = gs.getSize()[i] - 1;
            else
                sz_g[i] = gs.getSize()[i];
        }
 
        CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);
 
        for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
            spacing[i] = cd.getCellBox().getHigh(i);
    }
 
public:
 
    void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
    {
        for (size_t i = 0 ; i < Sys_eqs::nvar ; i++)
            s_pos[i] = sp[i];
    }
 
    void computeStag()
    {
        typedef typename Sys_eqs::b_grid::value_type prp_type;
 
        openfpm::vector<comb<Sys_eqs::dims>> c_prp[prp_type::max_prop];
 
        stag_set_position<Sys_eqs::dims,prp_type> ssp(c_prp);
 
        boost::mpl::for_each_ref< boost::mpl::range_c<int,0,prp_type::max_prop> >(ssp);
    }
 
    const Padding<Sys_eqs::dims> & getPadding()
    {
        return pd;
    }
 
    const g_map_type & getMap()
    {
        return g_map;
    }
 
    FDScheme(const Ghost<Sys_eqs::dims,long int> & stencil,
             const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
             const typename Sys_eqs::b_grid & b_g)
    :pd({0,0,0},{0,0,0}),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
    {
        Initialize(domain);
    }
 
    FDScheme(Padding<Sys_eqs::dims> & pd,
             const Ghost<Sys_eqs::dims,long int> & stencil,
             const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
             const typename Sys_eqs::b_grid & b_g)
    :pd(pd),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
        {
        Initialize(domain);
    }
 
    template<typename T> void impose(const T & op,
                                     typename Sys_eqs::stype num,
                                     long int id,
                                     const long int (& start)[Sys_eqs::dims],
                                     const long int (& stop)[Sys_eqs::dims],
                                     bool skip_first = false)
    {
    grid_key_dx<Sys_eqs::dims> start_k;
    grid_key_dx<Sys_eqs::dims> stop_k;
 
    bool increment = false;
    if (skip_first == true)
    {
            start_k = grid_key_dx<Sys_eqs::dims>(start);
            stop_k = grid_key_dx<Sys_eqs::dims>(start);
 
            auto it = g_map.getSubDomainIterator(start_k,stop_k);
 
            if (it.isNext() == true)
                    increment = true;
    }
 
        // add padding to start and stop
        start_k = grid_key_dx<Sys_eqs::dims>(start);
        stop_k = grid_key_dx<Sys_eqs::dims>(stop);
 
        auto it = g_map.getSubDomainIterator(start_k,stop_k);
 
        if (increment == true)
        {++it;}
 
        constant_b b(num);
 
        impose_git_gmap(op,b,id,it);
 
    }
 
    void new_b()
    {row_b = 0;}
 
    void new_A()
    {row = 0;}
 
    template<unsigned int prp, typename b_term, typename iterator>
    void impose_dit_b(b_term & b_t,
                      const iterator & it_d,
                      long int id = 0)
    {
        grid_b<b_term,prp> b(b_t);
 
        impose_dit_b(b,id,it_d);
    }
 
    template<typename T> void impose_dit(const T & op ,
                                     typename Sys_eqs::stype num,
                                     long int id ,
                                     grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
    {
        constant_b b(num);
 
        impose_dit(op,b,id,it_d);
    }
 
    template<typename T, typename bop, typename iterator> void impose_git_gmap(const T & op ,
                                     bop num,
                                     long int id ,
                                     iterator & it)
    {
        openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
 
        grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
 
        std::unordered_map<long int,float> cols;
 
        // iterate all the grid points
        while (it.isNext())
        {
            // get the position
            auto key = it.get();
 
            // Calculate the non-zero colums
            T::value(g_map,key,gs,spacing,cols,1.0);
 
            // indicate if the diagonal has been set
            bool is_diag = false;
 
            // create the triplet
            for ( auto it = cols.begin(); it != cols.end(); ++it )
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().col() = it->first;
                trpl.last().value() = it->second;
 
                if (trpl.last().row() == trpl.last().col())
                    is_diag = true;
 
//              std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
            }
 
            // If does not have a diagonal entry put it to zero
            if (is_diag == false)
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
                trpl.last().value() = 0.0;
            }
 
            b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
 
            cols.clear();
 
            // if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//          T::position(key,gs,s_pos);
#endif
 
            ++row;
            ++row_b;
            ++it;
        }
    }
 
    template<typename T, typename bop, typename iterator> void impose_git(const T & op ,
                                     bop num,
                                     long int id ,
                                     const iterator & it_d,
                                     bool skip_first = false)
    {
        openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
 
        auto start = it_d.getStart();
        auto stop = it_d.getStop();
 
        auto itg = g_map.getSubDomainIterator(start,stop);
 
        auto it = it_d;
 
        if (skip_first == true)
        {++it;}
 
        grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
 
        std::unordered_map<long int,float> cols;
 
        // iterate all the grid points
        while (it.isNext())
        {
            // get the position
            auto key = it.get();
            auto keyg = itg.get();
 
            // Calculate the non-zero colums
            T::value(g_map,keyg,gs,spacing,cols,1.0);
 
            // indicate if the diagonal has been set
            bool is_diag = false;
 
            // create the triplet
            for ( auto it = cols.begin(); it != cols.end(); ++it )
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
                trpl.last().col() = it->first;
                trpl.last().value() = it->second;
 
                if (trpl.last().row() == trpl.last().col())
                    is_diag = true;
 
//              std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
            }
 
            // If does not have a diagonal entry put it to zero
            if (is_diag == false)
            {
                trpl.add();
                trpl.last().row() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
                trpl.last().col() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
                trpl.last().value() = 0.0;
            }
 
            b(g_map.template get<0>(keyg)*Sys_eqs::nvar + id) = num.get(key);
 
            cols.clear();
 
            // if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//          T::position(key,gs,s_pos);
#endif
 
            ++row;
            ++row_b;
            ++it;
            ++itg;
        }
    }
 
    template<unsigned int prp, typename T, typename b_term, typename iterator> void impose_dit(const T & op ,
                                     b_term & b_t,
                                     const iterator & it_d,
                                     long int id = 0)
    {
        grid_b<b_term,prp> b(b_t);
 
        impose_dit(op,b,id,it_d);
    }
 
    typename Sys_eqs::SparseMatrix_type A;
 
    typename Sys_eqs::SparseMatrix_type & getA()
    {
#ifdef SE_CLASS1
        consistency();
#endif
        A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar,
                 g_map.getLocalDomainSize()*Sys_eqs::nvar,
                 g_map.getLocalDomainSize()*Sys_eqs::nvar);
 
        return A;
 
    }
 
    typename Sys_eqs::Vector_type & getB()
    {
#ifdef SE_CLASS1
        consistency();
#endif
 
        return b;
    }
 
    template<unsigned int ... pos, typename Vct, typename Grid_dst>
    void copy(Vct & v,const long int (& start)[Sys_eqs_typ::dims], const long int (& stop)[Sys_eqs_typ::dims], Grid_dst & g_dst)
    {
        if (is_grid_staggered<Sys_eqs>::value())
        {
            if (g_dst.is_staggered() == true)
            {
                grid_key_dx<Grid_dst::dims> sr(start);
                grid_key_dx<Grid_dst::dims> st(stop);
                copy_staggered<Vct,Grid_dst,pos...>(v,g_dst,sr,st);
            }
            else
            {
                // Create a temporal staggered grid and copy the data there
                auto & g_map = this->getMap();
 
                // Convert the ghost in grid units
 
                Ghost<Grid_dst::dims,long int> g_int;
                for (size_t i = 0 ; i < Grid_dst::dims ; i++)
                {
                    g_int.setLow(i,g_map.getDecomposition().getGhost().getLow(i) / g_map.spacing(i));
                    g_int.setHigh(i,g_map.getDecomposition().getGhost().getHigh(i) / g_map.spacing(i));
                }
 
                staggered_grid_dist<Grid_dst::dims,typename Grid_dst::stype,typename Grid_dst::value_type,typename Grid_dst::decomposition::extended_type, typename Grid_dst::memory_type, typename Grid_dst::device_grid_type> stg(g_dst,g_int,this->getPadding());
                stg.setDefaultStagPosition();
 
                grid_key_dx<Grid_dst::dims> sr;
                grid_key_dx<Grid_dst::dims> st;
 
                for (int i = 0 ; i < Grid_dst::dims ; i++)
                {
                    sr.set_d(i,-1);
                    st.set_d(i,stop[i]);
                }
 
                copy_staggered<Vct,decltype(stg),pos...>(v,stg,sr,st);
 
                // sync the ghost and interpolate to the normal grid
                stg.template ghost_get<pos...>();
                stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
            }
        }
        else
        {
            copy_normal<Vct,Grid_dst,pos...>(v,g_dst);
        }
    }
 
    template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v, Grid_dst & g_dst)
    {
        long int start[Sys_eqs_typ::dims];
        long int stop[Sys_eqs_typ::dims];
 
        for (size_t i = 0 ; i < Sys_eqs_typ::dims ; i++)
        {
            start[i] = 0;
            stop[i] = g_dst.size(i);
        }
 
        this->copy<pos...>(v,start,stop,g_dst);
    }
};
 
#define EQS_FIELDS 0
#define EQS_SPACE 1
 
 
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */