CartDecomposition.hpp

/*
 * CartDecomposition.hpp
 *
 *  Created on: Oct 07, 2015
 *      Author: Pietro Incardona, Antonio Leo
 */

#ifndef CARTDECOMPOSITION_HPP
#define CARTDECOMPOSITION_HPP

#include "config.h"
#include <cmath>
#include "VCluster/VCluster.hpp"
#include "Graph/CartesianGraphFactory.hpp"
#include "Decomposition.hpp"
#include "Vector/map_vector.hpp"
#include <vector>
#include <initializer_list>
#include "SubdomainGraphNodes.hpp"
#include "dec_optimizer.hpp"
#include "Space/Shape/Box.hpp"
#include "Space/Shape/Point.hpp"
#include "NN/CellList/CellDecomposer.hpp"
#include <unordered_map>
#include "NN/CellList/CellList.hpp"
#include "Space/Ghost.hpp"
#include "common.hpp"
#include "ie_loc_ghost.hpp"
#include "ie_ghost.hpp"
#include "nn_processor.hpp"
#include "GraphMLWriter/GraphMLWriter.hpp"
#include "Distribution/ParMetisDistribution.hpp"
#include "Distribution/DistParMetisDistribution.hpp"
#include "Distribution/MetisDistribution.hpp"
#include "DLB/DLB.hpp"
#include "util/se_util.hpp"
#include "util/mathutil.hpp"
#include "CartDecomposition_ext.hpp"
#include "data_type/aggregate.hpp"
#include "Domain_NN_calculator_cart.hpp"
#include "cuda/CartDecomposition_gpu.cuh"
#include "Domain_icells_cart.hpp"

#define CARTDEC_ERROR 2000lu

/*! \brief It spread the sub-sub-domain on a regular cartesian grid of size dim
 *
 * \warning this function only guarantee that the division on each direction is
 *          2^n with some n and does not guarantee that the number of
 *          sub-sub-domain is preserved
 *
 * \param div number of division on each direction as output
 * \param n_sub number of sub-domain
 * \param dim_r dimension reduction
 *
 */
template<unsigned int dim> static void nsub_to_div2(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
	for (size_t i = 0; i < dim; i++)
	{
		if (i < dim_r)
		{div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim_r));}
		else
		{div[i] = 1;}
	}
}

/*! \brief It spread the sub-sub-domain on a regular cartesian grid of size dim
 *
 * \warning this function only guarantee that the division on each direction is
 *          2^n with some n and does not guarantee that the number of
 *          sub-sub-domain is preserved
 *
 * \param div number of division on each direction as output
 * \param n_sub number of sub-domain
 * \param dim_r dimension reduction
 *
 */
template<unsigned int dim> static void nsub_to_div(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
	for (size_t i = 0; i < dim; i++)
	{
		if (i < dim_r)
		{div[i] = std::floor(pow(n_sub, 1.0 / dim_r));}
		else
		{div[i] = 1;}
	}
}

#define COMPUTE_SKIN_SUB 1

/**
 * \brief This class decompose a space into sub-sub-domains and distribute them across processors
 *
 * \tparam dim is the dimensionality of the physical domain we are going to decompose.
 * \tparam T type of the space we decompose, Real, Integer, Complex ...
 * \tparam Memory Memory factory used to allocate memory
 * \tparam Distribution type of distribution, can be ParMetisDistribution or MetisDistribution
 *
 * Given an N-dimensional space, this class decompose the space into a Cartesian grid of small
 * sub-sub-domain. To each sub-sub-domain is assigned an id that identify at which processor is
 * assigned (in general the union of all the sub-sub-domain assigned to a processor is
 * simply connected space), a second step merge several sub-sub-domain with same id into bigger region
 *  sub-domain. Each sub-domain has an extended space called ghost part
 *
 * Assuming that VCluster.getProcessUnitID(), equivalent to the MPI processor rank, return the processor local
 * processor id, we define
 *
 * * local processor: processor rank
 * * local sub-domain: sub-domain given to the local processor
 * * external ghost box: (or ghost box) are the boxes that compose the ghost space of the processor, or the
 *   boxes produced expanding every local sub-domain by the ghost extension and intersecting with the sub-domain
 *   of the other processors
 * * Near processors are the processors adjacent to the local processor, where with adjacent we mean all the processor
 *   that has a non-zero intersection with the ghost part of the local processor, or all the processors that
 *   produce non-zero external boxes with the local processor, or all the processor that should communicate
 *   in case of ghost data synchronization
 * * internal ghost box: is the part of ghost of the near processor that intersect the space of the
 *       processor, or the boxes produced expanding the sub-domain of the near processors with the local sub-domain
 * * Near processor sub-domain: is a sub-domain that live in the a near (or contiguous) processor
 * * Near processor list: the list of all the near processor of the local processor (each processor has a list
 *                        of the near processor)
 * * Local ghosts internal or external are all the ghosts that does not involve inter-processor communications
 *
 * \see calculateGhostBoxes() for a visualization of internal and external ghost boxes
 *
 * ### Create a Cartesian decomposition object on a Box space, distribute, calculate internal and external ghost boxes
 *
 * \snippet CartDecomposition_unit_test.hpp Create CartDecomposition
 *
 */

template<unsigned int dim, typename T, typename Memory, template <typename> class layout_base, typename Distribution>
class CartDecomposition: public ie_loc_ghost<dim, T,layout_base, Memory>,
						 public nn_prcs<dim, T,layout_base,Memory>,
						 public ie_ghost<dim,T,Memory,layout_base>,
						 public domain_nn_calculator_cart<dim>,
						 public domain_icell_calculator<dim,T,layout_base,Memory>
{
public:

	//! Type of the domain we are going to decompose
	typedef T domain_type;

	//! It simplify to access the SpaceBox element
	typedef SpaceBox<dim, T> Box;

	//! This class is base of itself
	typedef CartDecomposition<dim,T,Memory,layout_base,Distribution> base_type;

	//! This class admit a class defined on an extended domain
	typedef CartDecomposition_ext<dim,T,Memory,layout_base,Distribution> extended_type;

protected:

	//! bool that indicate whenever the buffer has been already transfer to device
	bool host_dev_transfer = false;

	//! Indicate the communication weight has been set
	bool commCostSet = false;

	//! This is the key type to access  data_s, for example in the case of vector
	//! acc_key is size_t
	typedef typename openfpm::vector<SpaceBox<dim, T>,
			Memory,
			typename memory_traits_lin<SpaceBox<dim, T>>::type,
			memory_traits_lin,
			openfpm::vector_grow_policy_default,
			openfpm::vect_isel<SpaceBox<dim, T>>::value>::access_key acc_key;

	//! the set of all local sub-domain as vector
	openfpm::vector<SpaceBox<dim, T>,Memory,typename layout_base<SpaceBox<dim, T>>::type,layout_base> sub_domains;

	//! the remote set of all sub-domains as vector of 'sub_domains' vectors
	mutable openfpm::vector<Box_map<dim, T>,Memory,typename layout_base<Box_map<dim, T>>::type,layout_base> sub_domains_global;

	//! for each sub-domain, contain the list of the neighborhood processors
	openfpm::vector<openfpm::vector<long unsigned int> > box_nn_processor;

	//! Structure that contain for each sub-sub-domain box the processor id
	//! exist for efficient global communication
	CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> fine_s;

	//! Structure that store the cartesian grid information
	grid_sm<dim, void> gr;

	//! Structure that store the cartesian grid information
	grid_sm<dim, void> gr_dist;

	//! Structure that decompose the space into cells without creating them
	//! useful to convert positions to CellId or sub-domain id in this case
	CellDecomposer_sm<dim, T, shift<dim,T>> cd;

	//! rectangular domain to decompose
	::Box<dim,T> domain;

	//! Box Spacing
	T spacing[dim];

	//! Magnification factor between distribution and