Commit 6e8741cf authored by Pietro Incardona's avatar Pietro Incardona

Several changes for periodic boundary conditions

parent 7f66938e
#! /bin/bash
# check if the directory $1/EIGEN exist
if [ -d "$1/EIGEN" ]; then
echo "EIGEN already installed"
exit 0
fi
wget http://ppmcore.mpi-cbg.de/upload/eigen-3.2.7.tar.bz2
rm -rf eigen-eigen-b30b87236a1b
tar -xf eigen-3.2.7.tar.bz2
cd eigen-eigen-b30b87236a1b
mkdir $1/EIGEN/
mv Eigen $1/EIGEN/Eigen
cd ..
rm -rf eigen-eigen-b30b87236a1b
#! /bin/bash
# check if the directory $1/OPENBLAS exist
if [ -d "$1/OPENBLAS" ]; then
echo "OPENBLAS already installed"
exit 0
fi
wget http://ppmcore.mpi-cbg.de/upload/OpenBLAS-0.2.15.tar.gz
rm -rf OpenBLAS-0.2.15
tar -xf OpenBLAS-0.2.15.tar.gz
cd OpenBLAS-0.2.15
# configuration
make
mkdir $1/OPENBLAS
make install PREFIX=$1/OPENBLAS
rm -rf OpenBLAS-0.2.15
#! /bin/bash
# check if the directory $1/SUITESPARSE exist
if [ -d "$1/SUITESPARSE" ]; then
echo "SUITESPARSE already installed"
exit 0
fi
#wget http://ppmcore.mpi-cbg.de/upload/SuiteSparse-4.4.5.tar.gz
rm -rf SuiteSparse
tar -xf SuiteSparse-4.4.5.tar.gz
cd SuiteSparse
# configuration
sed -i "/INSTALL_LIB\s=\s\/usr\/local\/lib/c\INSTALL_LIB = $1\/SUITESPARSE\/lib" SuiteSparse_config/SuiteSparse_config.mk
sed -i "/INSTALL_INCLUDE\s=\s\/usr\/local\/include/c\INSTALL_INCLUDE = $1\/SUITESPARSE\/include" SuiteSparse_config/SuiteSparse_config.mk
sed -i "/\sLAPACK\s=\s-llapack/c\LAPACK = " SuiteSparse_config/SuiteSparse_config.mk
sed -i "/\sBLAS\s=\s\-lopenblas/c\BLAS = -L$1/OPENBLAS/lib -lopenblas" SuiteSparse_config/SuiteSparse_config.mk
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$1/OPENBLAS/lib"
make
mkdir $1/SUITESPARSE
mkdir $1/SUITESPARSE/lib
mkdir $1/SUITESPARSE/include
make install
rm -rf SuiteSparse
......@@ -29,6 +29,7 @@
#include "ie_ghost.hpp"
#include "nn_processor.hpp"
#include "util/se_util.hpp"
#include "util/mathutil.hpp"
#define CARTDEC_ERROR 2000lu
......@@ -312,6 +313,9 @@ private:
// Initialize the geo_cell structure
ie_ghost<dim,T>::Initialize_geo_cell(domain,div,orig);
// Initialize shift vectors
ie_ghost<dim,T>::generateShiftVectors(domain);
}
/*! \brief Create the subspaces that decompose your domain
......@@ -465,6 +469,26 @@ public:
}
};
/*! \brief class to select the returned id by ghost_processorID
*
*/
class shift_id
{
public:
/*! \brief Return the shift id
*
* \param p structure containing the id informations
* \param b_id box_id
*
* \return shift_id id
*
*/
inline static size_t id(p_box<dim,T> & p, size_t b_id)
{
return p.shift_id;
}
};
/*! It calculate the internal ghost boxes
*
* Example: Processor 10 calculate
......@@ -581,9 +605,7 @@ p1[0]<-----+ +----> p2[0]
ie_ghost<dim,T>::create_box_nn_processor_int(v_cl,ghost,sub_domains,box_nn_processor,*this);
// ebox must come after ibox (in this case)
ie_loc_ghost<dim,T>::create_loc_ghost_ibox(ghost,sub_domains);
ie_loc_ghost<dim,T>::create_loc_ghost_ebox(ghost,sub_domains);
ie_loc_ghost<dim,T>::create(sub_domains,domain,ghost,bc);
// get the smallest sub-domain dimension on each direction
for (size_t i = 0 ; i < dim ; i++)
......@@ -622,11 +644,11 @@ p1[0]<-----+ +----> p2[0]
cart.ss_box = ss_box;
cart.ghost = g;
nn_prcs<dim,T>::create(box_nn_processor, sub_domains);
nn_prcs<dim,T>::applyBC(domain,ghost,bc);
(static_cast<nn_prcs<dim,T> &>(cart)).create(box_nn_processor, sub_domains);
(static_cast<nn_prcs<dim,T> &>(cart)).applyBC(domain,ghost,bc);
calculateGhostBoxes();
Initialize_geo_cell_lists();
cart.Initialize_geo_cell_lists();
cart.calculateGhostBoxes();
return cart;
}
......@@ -747,14 +769,14 @@ p1[0]<-----+ +----> p2[0]
return fine_s.get(cd.getCell(p));
}
/*! \brief Get the smallest subdivision of the domain on each direction
/*! \brief Given a point return in which processor the particle should go
*
* \return a box p1 is set to zero
* \return processorID
*
*/
const ::Box<dim,T> & getSmallestSubdivision()
size_t inline processorID(const Point<dim,T> &p) const
{
return ss_box;
return fine_s.get(cd.getCell(p));
}
/*! \brief Given a point return in which processor the particle should go
......@@ -762,19 +784,28 @@ p1[0]<-----+ +----> p2[0]
* \return processorID
*
*/
size_t inline processorID(const T (&p)[dim]) const
{
return fine_s.get(cd.getCell(p));
}
/*! \brief Get the smallest subdivision of the domain on each direction
*
* \return a box p1 is set to zero
*
*/
const ::Box<dim,T> & getSmallestSubdivision()
{
return ss_box;
}
/*! \brief Set the parameter of the decomposition
*
* \param div_ storing into how many domain to decompose on each dimension
* \param domain_ domain to decompose
*
*/
void setParameters(const size_t (& div_)[dim], Domain<dim,T> domain_, const size_t (& bc)[dim] ,Ghost<dim,T> ghost = Ghost<dim,T>())
void setParameters(const size_t (& div_)[dim], Domain<dim,T> domain_, const size_t (& bc)[dim] ,const Ghost<dim,T> & ghost)
{
// set the boundary conditions
for (size_t i = 0 ; i < dim ; i++)
......@@ -854,6 +885,8 @@ p1[0]<-----+ +----> p2[0]
}
/*! \brief Check if the particle is local
*
* \warning if the particle id outside the domain the result is unreliable
*
* \param p object position
*
......@@ -866,6 +899,8 @@ p1[0]<-----+ +----> p2[0]
}
/*! \brief Check if the particle is local
*
* \warning if the particle id outside the domain the result is unreliable
*
* \param p object position
*
......@@ -877,6 +912,56 @@ p1[0]<-----+ +----> p2[0]
return processorID(pos) == v_cl.getProcessUnitID();
}
/*! \brief Check if the particle is local considering boundary conditions
*
* \warning if the particle id outside the domain and non periodic the result
* is unreliable
*
*
* \param p object position
*
* \return true if it is local
*
*/
template<typename Mem> bool isLocalBC(const encapc<1, Point<dim,T>, Mem> p, const size_t (& bc)[dim]) const
{
Point<dim,T> pt = p;
for (size_t i = 0 ; i < dim ; i++)
{
if (bc[i] == PERIODIC)
pt.get(i) = openfpm::math::periodic_l(p.template get<0>()[i],domain.getHigh(i),domain.getLow(i));
}
return processorID<Mem>(pt) == v_cl.getProcessUnitID();
}
/*! \brief Check if the particle is local considering boundary conditions
*
* \param p object position
*
* \return true if it is local
*
*/
bool isLocalBC(const T (&p)[dim], const size_t (& bc)[dim]) const
{
Point<dim,T> pt = p;
for (size_t i = 0 ; i < dim ; i++)
{
if (pt.get(0) < -0.1)
{
int debug = 0;
debug++;
}
if (bc[i] == PERIODIC)
pt.get(i) = openfpm::math::periodic_l(p[i],domain.getHigh(i),domain.getLow(i));
}
return processorID(pt) == v_cl.getProcessUnitID();
}
/*! \brief Return the bounding box containing union of all the sub-domains for the local processor
*
* \return The bounding box
......
......@@ -175,10 +175,12 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_periodic_test)
}
// Check the consistency
bool val = dec.check_consistency();
BOOST_REQUIRE_EQUAL(val,true);
// write the decomposition
dec.write("test_out");
// We duplicate the decomposition
CartDecomposition<3,float> dec2 = dec.duplicate();
dec2.check_consistency();
......@@ -204,7 +206,6 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_periodic_test)
// Check that g3 is equal to dec2 with the exception of the ghost part
ret = dec3.is_equal_ng(dec2);
BOOST_REQUIRE_EQUAL(ret,true);
}
BOOST_AUTO_TEST_SUITE_END()
......
......@@ -13,26 +13,69 @@
#include "Vector/map_vector.hpp"
/*! \brief for each sub-domain box sub contain the real the sub-domain id
*
* When we apply boundary conditions real sub-domains are copied along the border
* sub, contain the id of the real sub_domain
*
* \tparam dim Dimensionality of the box
* \tparam T in witch space this box live
*
*/
template<unsigned int dim, typename T>
struct Box_loc_sub
{
Box<dim,T> bx;
// Domain id
size_t sub;
// in witch sector this sub-domain live, when
comb<dim> cmb;
Box_loc_sub()
{
cmb.zero();
};
Box_loc_sub(const Box<dim,T> & bx, size_t sub, const comb<dim> & cmb)
:bx(bx),sub(sub),cmb(cmb)
{};
template <typename Memory> Box_loc_sub(const Box_loc_sub<dim,T> & bls)
{
bx = bls.bx;
this->sub = bls.sub;
};
Box_loc_sub operator=(const Box<dim,T> & box)
{
::Box<dim,T>::operator=(box);
return *this;
}
};
/*! It contain a box definition and from witch sub-domain it come from (in the local processor)
* and an unique across adjacent processors (for communication)
*
* If the box come from the intersection of an expanded sub-domain and a sub-domain
*
* Assuming we are considering the adjacent processor i (0 to getNNProcessors())
* Assuming we are considering the near processors i (0 to getNNProcessors())
*
* ### external ghost box
*
* id = id_exp + N_non_exp + id_non_exp
* id = id_exp * N_non_exp + id_non_exp
*
* id_exp = the id in the vector proc_adj_box.get(i) of the expanded sub-domain
* id_exp = the id in the vector proc_adj_box.get(i) of the expanded sub-domain (sent local sub-domains)
*
* id_non_exp = the id in the vector nn_processor_subdomains[i] of the sub-domain
* id_non_exp = the id in the vector nn_processor_subdomains[i] of the sub-domain (received sub-domains from near processors)
*
* ### internal ghost box
*
* id = id_exp + N_non_exp + id_non_exp
* id = id_exp * N_non_exp + id_non_exp
*
* id_exp = the id in the vector nn_processor_subdomains[i] of the expanded sub-domain
*
......@@ -40,8 +83,10 @@
*
*/
template<unsigned int dim, typename T>
struct Box_sub : public Box<dim,T>
struct Box_sub
{
Box<dim,T> bx;
// Domain id
size_t sub;
......@@ -50,7 +95,7 @@ struct Box_sub : public Box<dim,T>
Box_sub operator=(const Box<dim,T> & box)
{
::Box<dim,T>::operator=(box);
bx = box;
return *this;
}
......@@ -60,17 +105,27 @@ struct Box_sub : public Box<dim,T>
//! Particular case for local internal ghost boxes
template<unsigned int dim, typename T>
struct Box_sub_k : public Box<dim,T>
struct Box_sub_k
{
Box<dim,T> bx;
// Domain id
size_t sub;
// Where this sub_domain live
comb<dim> cmb;
//! k \see getLocalGhostIBoxE
long int k;
Box_sub_k()
{
cmb.zero();
}
Box_sub_k operator=(const Box<dim,T> & box)
{
::Box<dim,T>::operator=(box);
bx = box;
return *this;
}
......@@ -129,9 +184,12 @@ struct N_box
// id of the processor in the nn_processor list (local processor id)
size_t id;
// adjacent processor sub-domains
// near processor sub-domains
typename openfpm::vector<::Box<dim,T>> bx;
// near processor sector position (or where they live outside the domain)
openfpm::vector<comb<dim>> pos;
//! Default constructor
N_box() {};
......@@ -208,9 +266,12 @@ struct p_box
::Box<dim,T> box;
//! local processor id
size_t lc_proc;
//! processor id
//! processor rank
size_t proc;
//! shift vector id
size_t shift_id;
/*! \brief Check if two p_box are the same
*
* \param pb box to check
......
/*
* decomposition_util_test.hpp
*
* Created on: Dec 16, 2015
* Author: i-bird
*/
#ifndef SRC_DECOMPOSITION_DECOMPOSITION_UTIL_TEST_HPP_
#define SRC_DECOMPOSITION_DECOMPOSITION_UTIL_TEST_HPP_
#include "VCluster.hpp"
openfpm::vector<SpaceBox<3,float>> create3Ddecomposition(Vcluster & vcl)
{
vb3.add(Box<3,float>({0.2,0.2,0.5},{1.0,0.5,1.0}));
vb3.add(Box<3,float>({0.0,0.0,0.5},{0.2,0.2,1.0}));
vb3.add(Box<3,float>({0.2,0.0,0.5},{0.5,0.2,1.0}));
vb3.add(Box<3,float>({0.5,0.0,0.5},{1.0,0.2,1.0}));
vb3.add(Box<3,float>({0.0,0.2,0.5},{0.2,0.5,1.0}));
vb3.add(Box<3,float>({0.0,0.5,0.5},{1.0,1.0,1.0}));
}
#endif /* SRC_DECOMPOSITION_DECOMPOSITION_UTIL_TEST_HPP_ */
This diff is collapsed.
......@@ -13,6 +13,7 @@
#include "Space/SpaceBox.hpp"
#include "common.hpp"
#include "VTKWriter.hpp"
#include "nn_processor.hpp"
/*! \brief structure that store and compute the internal and external local ghost box
*
......@@ -27,10 +28,8 @@ class ie_loc_ghost
{
openfpm::vector<lBox_dom<dim,T>> loc_ghost_box;
// Save the ghost boundaries
// Ghost<dim,T> ghost;
protected:
//! temporal added sub-domains
openfpm::vector<Box_loc_sub<dim,T>> sub_domains_tmp;
/*! \brief Create the external local ghost boxes
*
......@@ -38,11 +37,8 @@ protected:
* \param local sub-domain
*
*/
void create_loc_ghost_ebox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains)
void create_loc_ghost_ebox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains, openfpm::vector<Box_loc_sub<dim,T>> & sub_domains_prc)
{
// Save the ghost
// this->ghost = ghost;
loc_ghost_box.resize(sub_domains.size());
// For each sub-domain
......@@ -54,30 +50,29 @@ protected:
sub_with_ghost.enlarge(ghost);
// intersect with the other local sub-domains
for (size_t j = 0 ; j < sub_domains.size() ; j++)
for (size_t j = 0 ; j < sub_domains_prc.size() ; j++)
{
if (i == j)
continue;
size_t rj = sub_domains_prc.get(j).sub;
::Box<dim,T> bi;
bool intersect = sub_with_ghost.Intersect(::SpaceBox<dim,T>(sub_domains.get(j)),bi);
bool intersect = sub_with_ghost.Intersect(sub_domains_prc.get(j).bx,bi);
if (intersect == true)
{
Box_sub<dim,T> b;
b.sub = j;
b.sub = rj;
b = bi;
// local external ghost box
loc_ghost_box.get(i).ebx.add(b);
// search this box in the internal box of the sub-domain j
for (size_t k = 0; k < loc_ghost_box.get(j).ibx.size() ; k++)
for (size_t k = 0; k < loc_ghost_box.get(rj).ibx.size() ; k++)
{
if (loc_ghost_box.get(j).ibx.get(k).sub == i)
if (loc_ghost_box.get(rj).ibx.get(k).sub == i && loc_ghost_box.get(rj).ibx.get(k).cmb == sub_domains_prc.get(j).cmb.operator-())
{
loc_ghost_box.get(j).ibx.get(k).k = loc_ghost_box.get(i).ebx.size()-1;
loc_ghost_box.get(rj).ibx.get(k).k = loc_ghost_box.get(i).ebx.size()-1;
break;
}
}
......@@ -91,7 +86,7 @@ protected:
* \param ghost margin to enlarge
*
*/
void create_loc_ghost_ibox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains)
void create_loc_ghost_ibox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains, openfpm::vector<Box_loc_sub<dim,T>> & sub_domains_prc)
{
loc_ghost_box.resize(sub_domains.size());
......@@ -99,12 +94,10 @@ protected:
for (size_t i = 0 ; i < sub_domains.size() ; i++)
{
// intersect with the others local sub-domains
for (size_t j = 0 ; j < sub_domains.size() ; j++)
for (size_t j = 0 ; j < sub_domains_prc.size() ; j++)
{
if (i == j)
continue;
SpaceBox<dim,T> sub_with_ghost = sub_domains.get(j);
SpaceBox<dim,T> sub_with_ghost = sub_domains_prc.get(j).bx;
size_t rj = sub_domains_prc.get(j).sub;
// enlarge the sub-domain with the ghost
sub_with_ghost.enlarge(ghost);
......@@ -115,9 +108,10 @@ protected:
if (intersect == true)
{
Box_sub_k<dim,T> b;
b.sub = j;
b.sub = rj;
b = bi;
b.k = -1;
b.cmb = sub_domains_prc.get(j).cmb;
loc_ghost_box.get(i).ibx.add(b);
}
......@@ -125,8 +119,132 @@ protected:
}
}
/*! \brief In case of periodic boundary conditions we have to add boxes
* at the borders
*
* \param list of sub-domains
* \param domain Domain box
* \param boundary boundary conditions
* \param ghost ghost part
*
*/
void applyBC(openfpm::vector<Box_loc_sub<dim,T>> & sub_domains, const Box<dim,T> & domain, const Ghost<dim,T> & ghost, const size_t (&bc)[dim])
{
HyperCube<dim> hyp;
// first we create boxes at the border of the domain used to detect the sub-domain
// that must be adjusted, each of this boxes define a shift in case of periodic boundary condition
for (long int i = dim-1 ; i >= 0 ; i--)
{
std::vector<comb<dim>> cmbs = hyp.getCombinations_R(i);
for (size_t j = 0 ; j < cmbs.size() ; j++)
{
if (nn_prcs<dim,T>::check_valid(cmbs[j],bc) == false)
continue;
Box<dim,T> bp;
Point<dim,T> shift;
for (size_t k = 0 ; k < dim ; k++)
{
switch (cmbs[j][k])
{
case 1:
bp.setLow(k,domain.getHigh(k)+ghost.getLow(k));
bp.setHigh(k,domain.getHigh(k));
shift.get(k) = -domain.getHigh(k);
break;
case 0:
bp.setLow(k,domain.getLow(k));
bp.setHigh(k,domain.getHigh(k));
shift.get(k) = 0;
break;
case -1:
bp.setLow(k,domain.getLow(k));
bp.setHigh(k,ghost.getHigh(k));
shift.get(k) = domain.getHigh(k);
break;
}
}
// Detect all the sub-domain involved, shift them and add to the list
// Detection is performed intersecting the sub-domains with the ghost
// parts near the domain borders
for (size_t k = 0 ; k < sub_domains.size() ; k++)
{
Box<dim,T> sub = sub_domains.get(k).bx;
Box<dim,T> b_int;
if (sub.Intersect(bp,b_int) == true)