SpaceDistribution.hpp

/*
 * SpaceDistribution.hpp
 *
 *  Created on: Dec 3, 2016
 *      Author: i-bird
 */
 
#ifndef SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_
#define SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_
 
#include "util/mathutil.hpp"
#include "NN/CellList/CellDecomposer.hpp"
#include "Grid/grid_key_dx_iterator_hilbert.hpp"
 
template<unsigned int dim, typename T>
class SpaceDistribution
{
    Vcluster<> & v_cl;
 
    grid_sm<dim, void> gr;
 
    Box<dim, T> domain;
 
    Graph_CSR<nm_v<dim>, nm_e> gp;
 
 
public:
 
    SpaceDistribution(Vcluster<> & v_cl)
    :v_cl(v_cl)
    {
    }
 
    SpaceDistribution(const ParMetisDistribution<dim,T> & pm)
    :v_cl(pm.v_cl)
    {
        this->operator=(pm);
    }
 
    SpaceDistribution(SpaceDistribution<dim,T> && pm)
    :v_cl(pm.v_cl)
    {
        this->operator=(pm);
    }
 
    void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
    {
        size_t bc[dim];
 
        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = NON_PERIODIC;
 
        // Set grid and domain
        gr = grid;
        domain = dom;
 
        // Create a cartesian grid graph
        CartesianGraphFactory<dim, Graph_CSR<nm_v<dim>, nm_e>> g_factory_part;
        gp = g_factory_part.template construct<NO_EDGE, nm_v_id, T, dim - 1, 0>(gr.getSize(), domain, bc);
 
        // Init to 0.0 axis z (to fix in graphFactory)
        if (dim < 3)
        {
            for (size_t i = 0; i < gp.getNVertex(); i++)
                gp.vertex(i).template get<nm_v_x>()[2] = 0.0;
        }
        for (size_t i = 0; i < gp.getNVertex(); i++)
            gp.vertex(i).template get<nm_v_global_id>() = i;
 
    }
 
    Graph_CSR<nm_v<dim>, nm_e> & getGraph()
    {
        return gp;
    }
 
    void decompose()
    {
        // Get the number of processing units
        size_t Np = v_cl.getProcessingUnits();
 
        // Calculate the best number of sub-domains for each
        // processor
        size_t N_tot = gr.size();
        size_t N_best_each = N_tot / Np;
        size_t N_rest = N_tot % Np;
 
        openfpm::vector<size_t> accu(Np);
        accu.get(0) = N_best_each + ((0 < N_rest)?1:0);
        for (size_t i = 1 ; i < Np ; i++)
            accu.get(i) = accu.get(i-1) + N_best_each + ((i < N_rest)?1:0);
 
        // Get the maximum along dimensions and take the smallest n number
        // such that 2^n < m. n it will be order of the hilbert curve
 
        size_t max = 0;
 
        for (size_t i = 0; i < dim ; i++)
        {
            if (max < gr.size(i))
                max = gr.size(i);
        }
 
        // Get the order of the hilbert-curve
        size_t order = openfpm::math::log2_64(max);
        if (1ul << order < max)
            order += 1;
 
        size_t n = 1 << order;
 
        // Create the CellDecomoser
 
        CellDecomposer_sm<dim,T> cd_sm;
        cd_sm.setDimensions(domain, gr.getSize(), 0);
 
        // create the hilbert curve
 
        //hilbert curve iterator
        grid_key_dx_iterator_hilbert<dim> h_it(order);
 
        T spacing[dim];
 
        // Calculate the hilbert curve spacing
        for (size_t i = 0 ; i < dim ; i++)
            spacing[i] = (domain.getHigh(i) - domain.getLow(i)) / n;
 
        // Small grid to detect already assigned sub-sub-domains
        grid_cpu<dim,aggregate<long int>> g(gr.getSize());
        g.setMemory();
 
        // Reset the grid to -1
        grid_key_dx_iterator<dim> it(gr);
        while (it.isNext())
        {
            auto key = it.get();
 
            g.template get<0>(key) = -1;
 
            ++it;
        }
 
        // Go along the hilbert-curve and divide the space
 
        size_t proc_d = 0;
        size_t ele_d = 0;
        while (h_it.isNext())
        {
          auto key = h_it.get();
 
          // Point p
          Point<dim,T> p;
 
          for (size_t i = 0 ; i < dim ; i++)
              p.get(i) = key.get(i) * spacing[i] + spacing[i] / 2;
 
          grid_key_dx<dim> sp = cd_sm.getCellGrid(p);
 
          if (g.template get<0>(sp) == -1)
          {
              g.template get<0>(sp) = proc_d;
              ele_d++;
 
              if (ele_d >= accu.get(proc_d))
                  proc_d++;
          }
 
          ++h_it;
        }
 
        // Fill from the grid to the graph
 
        // Reset the grid to -1
        grid_key_dx_iterator<dim> it2(gr);
        while (it2.isNext())
        {
            auto key = it2.get();
 
            gp.template vertex_p<nm_v_proc_id>(gr.LinId(key)) = g.template get<0>(key);
 
            ++it2;
        }
 
        return;
    }
 
    void refine()
    {
        std::cout << __FILE__ << ":" << __LINE__ << " You are trying to dynamicaly balance a fixed decomposition, this operation has no effect" << std::endl;
    }
 
    float getUnbalance()
    {
        return gr.size() % v_cl.getProcessingUnits();
    }
 
    void getSubSubDomainPosition(size_t id, T (&pos)[dim])
    {
#ifdef SE_CLASS1
        if (id >= gp.getNVertex())
            std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
#endif
 
        // Copy the geometrical informations inside the pos vector
        pos[0] = gp.vertex(id).template get<nm_v_x>()[0];
        pos[1] = gp.vertex(id).template get<nm_v_x>()[1];
        if (dim == 3)
            pos[2] = gp.vertex(id).template get<nm_v_x>()[2];
    }
 
    inline void setComputationCost(size_t id, size_t weight)
    {
        std::cout << __FILE__ << ":" << __LINE__ << " You are trying to set the computation cost on a fixed decomposition, this operation has no effect" << std::endl;
    }
 
    bool weightsAreUsed()
    {
        return false;
    }
 
    size_t getSubSubDomainComputationCost(size_t id)
    {
        return 1.0;
    }
 
    size_t getProcessorLoad()
    {
        // Get the number of processing units
        size_t Np = v_cl.getProcessingUnits();
 
        // Calculate the best number of sub-domains for each
        // processor
        size_t N_tot = gr.size();
        size_t N_best_each = N_tot / Np;
        size_t N_rest = N_tot % Np;
 
        if (v_cl.getProcessUnitID() < N_rest)
            N_best_each += 1;
 
        return N_best_each;
    }
 
    void setMigrationCost(size_t id, size_t migration)
    {
    }
 
    void setCommunicationCost(size_t v_id, size_t e, size_t communication)
    {
    }
 
    size_t getNSubSubDomains()
    {
        return gp.getNVertex();
    }
 
    size_t getNSubSubDomainNeighbors(size_t id)
    {
        return gp.getNChilds(id);
    }
 
    void write(const std::string & file)
    {
        VTKWriter<Graph_CSR<nm_v<dim>, nm_e>, VTK_GRAPH> gv2(gp);
        gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + file + ".vtk");
    }
 
    const SpaceDistribution<dim,T> & operator=(const SpaceDistribution<dim,T> & dist)
    {
        gr = dist.gr;
        domain = dist.domain;
        gp = dist.gp;
 
        return *this;
    }
 
    const SpaceDistribution<dim,T> & operator=(SpaceDistribution<dim,T> && dist)
    {
        //v_cl = dist.v_cl;
        gr = dist.gr;
        domain = dist.domain;
        gp.swap(dist.gp);
 
        return *this;
    }
 
    size_t get_ndec()
    {
        return 0;
    }
};
 
 
#endif /* SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_ */