diff --git a/CHANGELOG.md b/CHANGELOG.md index c04b0b9bca0253646874ca0acf04d5dd8b46ff20..54f0e06a995be1881ecfb4a5faa00bb3d713ed79 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,19 @@ # Change Log All notable changes to this project will be documented in this file. +## [0.2.0] - 2015-02-17 +### Added +- PSE 1D example with multiple precision +- Plot example for GoogleChart plotting +- Distributed data structure now support 128bit floating point precision (on Beta) +- OpenFPM support for Microsoft Windows (Cygwin) compilation + +### Fixed +- Detection 32 bit system and report as an error + +### Changed +- Nothing to report + ## [0.1.0] - 2015-02-05 ### Added - PSE 1D example diff --git a/Makefile.am b/Makefile.am index 2a675a02ee48361d4af1c61b67be361fa39e71af..a6fbc5374b3fd1f026c8d0ad28d5b5b393c69a7d 100755 --- a/Makefile.am +++ b/Makefile.am @@ -1,3 +1,3 @@ -SUBDIRS = src vtk openfpm_data openfpm_io openfpm_devices openfpm_vcluster +SUBDIRS = src vtk openfpm_data openfpm_io openfpm_devices openfpm_vcluster openfpm_numerics bin_PROGRAMS = diff --git a/configure.ac b/configure.ac index 9945dfa989743a65af10646fd488da7e8fc87812..295e144b2f6fab33ac60d939b2644bf5f737ba1a 100755 --- a/configure.ac +++ b/configure.ac @@ -40,7 +40,17 @@ m4_ifdef([AX_EIGEN],,[m4_include([m4/ax_eigen.m4])]) m4_ifdef([AX_LIB_HDF5],,[m4_include([m4/ax_lib_hdf5.m4])]]) m4_ifdef([AX_H5PART],,[m4_include([m4/ax_h5part.m4])]) -CXXFLAGS+=" --std=c++11 " +case $host_os in + *cygwin*) + # Do something specific for cygwin + CXXFLAGS+=" --std=gnu++11 " + ;; + *) + #Default Case + CXXFLAGS+=" --std=c++11 " + ;; +esac + NVCCFLAGS=" " INCLUDES_PATH=" " diff --git a/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp b/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp index 3fd936b4b71a1900c620fba3e2ed7deb6bc443d1..cd16f24ea75324682efd284c1c71134d53c9cba4 100644 --- a/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp +++ b/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp @@ -11,7 +11,7 @@ #include "Vector/vector_dist.hpp" #include "Decomposition/CartDecomposition.hpp" -#include "Kernels.hpp" +#include "PSE/Kernels.hpp" #include "data_type/aggregate.hpp" #include <cmath> @@ -63,10 +63,10 @@ int main(int argc, char* argv[]) // // Number of particles - const size_t Npart = 1000; + const size_t Npart = 125; // Number of padding particles (At least) - const size_t Npad = 20; + const size_t Npad = 40; // The domain Box<1,double> box({0.0},{4.0}); @@ -212,7 +212,7 @@ int main(int argc, char* argv[]) // get and construct the Cell list Ghost<1,double> gp(enlarge); - auto cl = vd.getCellList(8*eps,gp); + auto cl = vd.getCellList(12*eps,gp); // Maximum infinity norm double linf = 0.0; diff --git a/example/SE/0_classes/Makefile b/example/SE/0_classes/Makefile index 172ffbd665222f8f35a17e1b68d2d628a8f07208..3cd88b493e5e1adc65ccc44010a4aee0c9036fab 100644 --- a/example/SE/0_classes/Makefile +++ b/example/SE/0_classes/Makefile @@ -7,7 +7,7 @@ LDIR = OBJ = main.o %.o: %.cpp - $(CC) -O0 -g3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH) + $(CC) -O3 -g3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH) se_classes: $(OBJ) $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS) diff --git a/m4/ax_suitesparse.m4 b/m4/ax_suitesparse.m4 index f4ff1309ede91dd1ebee35dd1e8d19960d1f9969..ea32ce0ed11348b27b1d17511ceed7c7f564e849 100644 --- a/m4/ax_suitesparse.m4 +++ b/m4/ax_suitesparse.m4 @@ -87,10 +87,11 @@ esac # First, check SUITESPARSE_LIBS environment variable if test "x$SUITESPARSE_LIBS" != x; then - save_LIBS="$LIBS"; LIBS="$SUITESPARSE_LIBS -lumfpack -lsuitesparseconfig -lm $RT_LIB" + save_LIBS="$LIBS"; LIBS="$SUITESPARSE_LIBS -lumfpack -lamd -lbtf -lcamd -lccolamd -lcholmod -lcolamd -lcxsparse -lklu -ldl -lrbio -lspqr -lsuitesparseconfig -lm $RT_LIB" AC_MSG_CHECKING([for umf_l_malloc]) AC_TRY_LINK_FUNC(umf_l_malloc, [ax_suitesparse_ok=yes - SUITESPARSE_LIB="$SUITESPARSE_LIBS -lamd -lbtf -lcamd -lccolamd -lcholmod -lcolamd -lcxsparse -lklu -ldl -lrbio -lspqr -lsuitesparseconfig -lumfpack"], [SUITRSPARSE_LIBS=""]) + SUITESPARSE_LIBS="$SUITESPARSE_LIBS -lumfpack -lamd -lbtf -lcamd -lccolamd -lcholmod -lcolamd -lcxsparse -lklu -ldl -lrbio -lspqr -lsuitesparseconfig" + ], [SUITRSPARSE_LIBS=""]) AC_MSG_RESULT($ax_suitesparse_ok) LIBS="$save_LIBS" if test $ax_suitesparse_ok = no; then @@ -100,12 +101,18 @@ if test "x$SUITESPARSE_LIBS" != x; then CFLAGS=$SUITESPARSE_INCLUDE AC_CHECK_HEADER(umfpack.h,[],[SUITESPARSE_INCLUDE="" ax_suitesparse_ok=no]) + CFLAGS="$old_CFLAGS" else - AC_CHECK_LIB(umfpack,umf_l_alloc,[SUITESPARSE_LIB="$SUITESPARSE_LIB -lamd -lbtf -lcamd -lccolamd -lcholmod -lcolamd -lcxsparse -lklu -ldl -lrbio -lspqr -lsuitesparseconfig -lumfpack"],[SUITESPARSE_LIB=""]) + AC_CHECK_LIB(umfpack,umf_l_alloc,[SUITESPARSE_LIBS="$SUITESPARSE_LIBS -lumfpack -lamd -lbtf -lcamd -lccolamd -lcholmod -lcolamd -lcxsparse -lklu -ldl -lrbio -lspqr -lsuitesparseconfig"],[ + SUITESPARSE_LIBS="" + ax_suitesparse_ok=no + ]) old_CFLAGS="$CFLAGS" AC_CHECK_HEADER(umfpack.h,[],[SUITESPARSE_INCLUDE="" ax_suitesparse_ok=no]) + + CFLAGS="$old_CFLAGS" fi diff --git a/script/discover_os b/script/discover_os index 722665c7a1443e232862d748320d40962cfda035..4abdd2de3e9737ab0f71a87527572c74bfe07539 100755 --- a/script/discover_os +++ b/script/discover_os @@ -15,8 +15,6 @@ platform=unknown platform=osx elif [[ "$OSTYPE" == "cygwin" ]]; then echo -e "We are on\033[1;34m CYGWIN \033[0m" - echo "This platform is not supported" - exit 1 elif [[ "$OSTYPE" == "msys" ]]; then echo -e "We are on\033[1;34m Microsoft Window \033[0m" echo "This platform is not supported" diff --git a/script/install_MPI.sh b/script/install_MPI.sh index e48c7275321c74b2b65260d6f8594d89e0aa5ce6..bc829f3bb4c336082513a2ae9427661ac234ddf0 100755 --- a/script/install_MPI.sh +++ b/script/install_MPI.sh @@ -10,6 +10,20 @@ fi wget http://ppmcore.mpi-cbg.de/upload/openmpi-1.8.7.tar.bz2 tar -xvf openmpi-1.8.7.tar.bz2 cd openmpi-1.8.7 + +# +# --disable-mca-dso \ +# --disable-sysv-shmem \ +# --enable-cxx-exceptions \ +# --with-threads=posix \ +# --without-cs-fs \ +# --with-mpi-param_check=always \ +# --enable-contrib-no-build=vt,libompitrace \ +# +#--enable-mca-no-build=paffinity,installdirs-windows,timer-windows,shmem-sysv +# +# + sh ./configure --prefix=$1/MPI --enable-opal-multi-threads --enable-mpi-f90 $2 $3 make -j 4 mkdir $1/MPI diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp index 2751dedf41fa10cb04854f9f4936c4c14e1116f5..217cfc1e36a03961a4cbc6c915a09b1b720024e3 100644 --- a/src/Decomposition/CartDecomposition.hpp +++ b/src/Decomposition/CartDecomposition.hpp @@ -351,6 +351,122 @@ private: } } + + /*! \brief It copy the sub-domains into another CartesianDecomposition object extending them + * + * \see duplicate (in case of extended domain) + * + * \param cart Cartesian decomposition object + * \param box Extended domain + * + */ + void extend_subdomains(CartDecomposition<dim,T> & cart, const ::Box<dim,T> & ext_dom) const + { + // Box + typedef ::Box<dim,T> b; + + cart.bbox = ext_dom; + cart.ss_box = ext_dom; + + for (size_t i = 0 ; i < sub_domains.size() ; i++) + { + ::Box<dim,T> box; + + // Calculate the extended box + for (size_t j = 0 ; j < dim ; j++) + { + if (sub_domains.template get<b::p1>(i)[j] == domain.getLow(j)) + box.setLow(j,ext_dom.getLow(j)); + else + box.setLow(j,sub_domains.template get<b::p1>(i)[j]); + + if (sub_domains.template get<b::p2>(i)[j] == domain.getHigh(j)) + box.setHigh(j,ext_dom.getHigh(j)); + else + box.setHigh(j,sub_domains.template get<b::p2>(i)[j]); + } + + // add the subdomain + cart.sub_domains.add(box); + + // Calculate the bound box + cart.bbox.enclose(box); + + // Create the smallest box contained in all sub-domain + cart.ss_box.contained(box); + } + } + + /*! \brief Extend the fines for the new Cartesian decomposition + * + * \param new_fines extended fine_s + * \param old_fines old fine_s + * + */ + void extend_fines(CartDecomposition<dim,T> & cart) const + { + // Extension, first we calculate the extensions of the new domain compared + // to the old one in cell units (each cell unit is a sub-sub-domain) + ::Box<dim,size_t> ext; + // Extension of the new fines structure + ::Box<dim,size_t> n_fines_ext; + // Extension of the old fines structure + ::Box<dim,size_t> o_fines_ext; + + size_t sz_new[dim]; + size_t sz_old[dim]; + + for (size_t i = 0; i < dim ; i++) + { + size_t p1 = (domain.getLow(i) - this->domain.getLow(i)) / cd.getCellBox().getHigh(i) + 1; + size_t p2 = (domain.getLow(i) - this->domain.getLow(i)) / cd.getCellBox().getHigh(i) + 1; + + ext.setLow(i,p1); + ext.setHigh(i,p2); + sz_new[i] = p1+p2+cd.getGrid().size(i); + sz_old[i] = cd.getGrid().size(i); + } + + grid_sm<dim,void> info_new(sz_new); + grid_sm<dim,void> info_old(sz_old); + + // resize the new fines + cart.fine_s.resize(info_new.size()); + + // we create an iterator that iterate across the full new fines + grid_key_dx_iterator<dim> fines_t(info_new); + + while (fines_t.isNext()) + { + auto key = fines_t.get(); + + // new_fines is bigger than old_fines structure + // out of bound key must be adjusted + // The adjustment produce a natural extension + // a representation can be seen in the figure of + // CartDecomposition duplicate function with extended domains + + grid_key_dx<dim> key_old; + for (size_t i = 0 ; i < dim ; i++) + { + key_old.set_d(i,(long int)key.get(i) - ext.getLow(i)); + if (key_old.get(i) < 0) + key_old.set_d(i,0); + else if(key_old.get(i) >= (long int)info_old.size(i) ) + key_old.set_d(i,info_old.size(i)-1); + } + + cart.fine_s.get(info_new.LinId(key)) = fine_s.get(info_old.LinId(key_old)); + + ++fines_t; + } + + cart.gr.setDimensions(sz_new); + + // the new extended CellDecomposer must be consistent with the old cellDecomposer. + cart.cd.setDimensions(cd,ext); + } + public: static constexpr int dims = dim; @@ -699,6 +815,83 @@ p1[0]<-----+ +----> p2[0] return cart; } + /*! \brief It create another object that contain the same decomposition information but with different ghost boxes and an extended domain + * + * The domain extension is produced extending the boxes at the border like in figure + * + * \verbatim + * ++--------------^--------^----------^----------+ +| | | | | +| A | E | F | N | +| +-----------------------------------+----> +| | | | | | | +| A | A | | F | | | +| | | | | | | +| | | E +----------+ N | N | +<--------------+ | | | | +| | | | | | | +| | | | G | | | +| | | | +----------> +| B | B | +----------+ | | +| | +--------+ | M | M | +| | | | H | | | +| | | +-----+----+----------> +<--------------+ D | | | | +| | | | I | L | L | +| C | C | | | | | +| | | | | | | +| +-----------------------------------+ | +| | | | | +| C | D | I | L | ++--------------v--------v-----v---------------+ + + * + * \endverbatim + * + * \param g ghost + * \param domain extended domain (MUST be extended) + * + * \return a duplicated decomposition with different ghost boxes and an extended domain + * + */ + CartDecomposition<dim,T,Memory> duplicate(const Ghost<dim,T> & g, const ::Box<dim,T> & ext_domain) const + { + CartDecomposition<dim,T,Memory> cart(v_cl); + + cart.box_nn_processor = box_nn_processor; + + // Calculate new sub-domains for extended domain + extend_subdomains(cart,ext_domain); + + // Calculate fine_s structure for the extended domain + // update the cell decomposer and gr + extend_fines(cart); + + // Get the old sub-sub-domain grid extension + + cart.domain = ext_domain; + + // spacing does not change + std::copy(spacing,spacing+3,cart.spacing); + + //! Runtime virtual cluster + cart.v_cl = v_cl; + + cart.ghost = g; + + for (size_t i = 0 ; i < dim ; i++) + cart.bc[i] = bc[i]; + + (static_cast<nn_prcs<dim,T> &>(cart)).create(cart.box_nn_processor, cart.sub_domains); + (static_cast<nn_prcs<dim,T> &>(cart)).applyBC(ext_domain,ghost,bc); + + cart.Initialize_geo_cell_lists(); + cart.calculateGhostBoxes(); + + return cart; + } + /*! \brief It create another object that contain the same information and act in the same way * * \return a duplicated decomposition diff --git a/src/Decomposition/CartDecomposition_unit_test.hpp b/src/Decomposition/CartDecomposition_unit_test.hpp index 834cad13ed3852d97248217514601dade15f8ae9..6d880e13318e210ceaa421723d4181ce20b76bb3 100644 --- a/src/Decomposition/CartDecomposition_unit_test.hpp +++ b/src/Decomposition/CartDecomposition_unit_test.hpp @@ -173,6 +173,7 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_periodic_test) } // Check the consistency + bool val = dec.check_consistency(); BOOST_REQUIRE_EQUAL(val,true); @@ -203,6 +204,141 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_periodic_test) BOOST_REQUIRE_EQUAL(ret,true); } +BOOST_AUTO_TEST_CASE( CartDecomposition_extend_test) +{ + // Initialize the global VCluster + init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv); + + // Vcluster + Vcluster & vcl = *global_v_cluster; + + //! [Create CartDecomposition] + CartDecomposition<3,float> dec(vcl); + + // Physical domain + Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0}); + Box<3,float> bulk({0.05,0.05,0.05},{0.95,0.95,0.95}); + Box<3,float> bulk_g({0.01,0.01,0.01},{0.95,0.95,0.95}); + size_t div[3]; + + // Get the number of processor and calculate the number of sub-domain + // for each processor (SUB_UNIT_FACTOR=64) + size_t n_proc = vcl.getProcessingUnits(); + size_t n_sub = n_proc * SUB_UNIT_FACTOR; + + // Set the number of sub-sub-domains on each dimension (in a scalable way) + for (int i = 0 ; i < 3 ; i++) + {div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/3));} + + // Define ghost size + Ghost<3,float> g(0.01); + + // Boundary conditions + size_t bc[] = {PERIODIC,PERIODIC,PERIODIC}; + + // Decompose + dec.setParameters(div,box,bc,g); + + // Internally create the ghosts + dec.calculateGhostBoxes(); + + // Duplicate the decomposition extending the domain + Box<3,float> box_ext({-0.1,-0.1,-0.1},{1.1,1.1,1.1}); + CartDecomposition<3,float> dec2 = dec.duplicate(g,box_ext); + + // than we create a grid iterator that iterate across the + + // Create a grid based iterator on the smaller decomposition + size_t sz[3] = {128,128,128}; + grid_dist_id_iterator_dec<CartDecomposition<3,float>> git(dec,sz); + + // iterate across the grid points and we check that the results are consistent + while (git.isNext()) + { + auto key = git.get(); + + Point<3,float> p = key.toPoint(); + + if (bulk_g.isInside(p) == false) + { + ++git; + continue; + } + + for (size_t i = 0 ; i < 3 ; i++) + p.get(i) = key.get(i) * git.getSpacing(i); + + const openfpm::vector<size_t> & pr_pid = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p); + const openfpm::vector<size_t> & pr2_pid = dec2.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p); + + BOOST_REQUIRE(pr_pid == pr2_pid); + + const openfpm::vector<size_t> & pr_bid = dec.template ghost_processorID<CartDecomposition<3,float>::box_id>(p); + const openfpm::vector<size_t> & pr2_bid = dec2.template ghost_processorID<CartDecomposition<3,float>::box_id>(p); + + BOOST_REQUIRE(pr_bid == pr2_bid); + + const openfpm::vector<size_t> & pr_lid = dec.template ghost_processorID<CartDecomposition<3,float>::lc_processor_id>(p); + const openfpm::vector<size_t> & pr2_lid = dec2.template ghost_processorID<CartDecomposition<3,float>::lc_processor_id>(p); + + BOOST_REQUIRE(pr_lid == pr2_lid); + + size_t n_test = 0; + + // Given a point p, if the point is near the boundary + if (bulk.isInside(p) == false) + { + // we can move p in one direction and check that dec2 return the same result + for (size_t i = 0 ; i < 3 ; i++) + { + Point<3,float> q = p; + + q.get(i) = p.get(i) + 0.05; + + if (box.isInside(q) == false) + { + const openfpm::vector<size_t> & pr3_pid = dec2.template ghost_processorID<CartDecomposition<3,float>::processor_id>(q); + const openfpm::vector<size_t> & pr3_bid = dec2.template ghost_processorID<CartDecomposition<3,float>::box_id>(q); + const openfpm::vector<size_t> & pr3_lid = dec2.template ghost_processorID<CartDecomposition<3,float>::lc_processor_id>(q); + BOOST_REQUIRE(pr2_pid == pr3_pid); + BOOST_REQUIRE(pr2_bid == pr3_bid); + BOOST_REQUIRE(pr2_lid == pr3_lid); + + n_test++; + } + + q = p; + q.get(i) = p.get(i) - 0.05; + + if (box.isInside(q) == false) + { + const openfpm::vector<size_t> & pr3_pid = dec2.template ghost_processorID<CartDecomposition<3,float>::processor_id>(q); + const openfpm::vector<size_t> & pr3_bid = dec2.template ghost_processorID<CartDecomposition<3,float>::box_id>(q); + const openfpm::vector<size_t> & pr3_lid = dec2.template ghost_processorID<CartDecomposition<3,float>::lc_processor_id>(q); + + if (pr2_pid != pr3_pid) + { + const openfpm::vector<size_t> & pr3_pid = dec2.template ghost_processorID<CartDecomposition<3,float>::processor_id>(q); + + int debug = 0; + debug++; + } + + BOOST_REQUIRE(pr2_pid == pr3_pid); + BOOST_REQUIRE(pr2_bid == pr3_bid); + BOOST_REQUIRE(pr2_lid == pr3_lid); + + n_test++; + } + } + } + + BOOST_REQUIRE(n_test != 0); + + ++git; + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/src/Decomposition/ie_loc_ghost.hpp b/src/Decomposition/ie_loc_ghost.hpp index ec2d08b6d57870766c063dfc33439879245d08fb..3772b98871b41ce2733144220f3a70bad42ba3c1 100644 --- a/src/Decomposition/ie_loc_ghost.hpp +++ b/src/Decomposition/ie_loc_ghost.hpp @@ -133,12 +133,11 @@ class ie_loc_ghost } } - /*! \brief In case of periodic boundary conditions we have to add boxes - * at the borders + /*! \brief In case of periodic boundary conditions we replicate the sub-domains at the border * * \param list of sub-domains * \param domain Domain box - * \param boundary boundary conditions + * \param boundary conditions * \param ghost ghost part * */ @@ -167,7 +166,7 @@ class ie_loc_ghost case 1: bp.setLow(k,domain.getHigh(k)+ghost.getLow(k)); bp.setHigh(k,domain.getHigh(k)); - shift.get(k) = -domain.getHigh(k); + shift.get(k) = -domain.getHigh(k)+domain.getLow(k); break; case 0: bp.setLow(k,domain.getLow(k)); @@ -177,7 +176,7 @@ class ie_loc_ghost case -1: bp.setLow(k,domain.getLow(k)); bp.setHigh(k,ghost.getHigh(k)); - shift.get(k) = domain.getHigh(k); + shift.get(k) = domain.getHigh(k)-domain.getLow(k); break; } } @@ -527,7 +526,10 @@ public: for (size_t j = 0 ; j < loc_ghost_box.get(i).ibx.size() ; j++) { if (loc_ghost_box.get(i).ibx.get(j).k == -1) + { + std::cout << "No ibx link" << "\n"; return false; + } } } @@ -536,9 +538,15 @@ public: for (size_t i = 0 ; i < loc_ghost_box.size() ; i++) { if (loc_ghost_box.get(i).ibx.size() == 0) + { + std::cout << "Zero ibx" << "\n"; return false; + } if (loc_ghost_box.get(i).ebx.size() == 0) + { + std::cout << "Zero ebx" << "\n"; return false; + } } } diff --git a/src/Decomposition/nn_processor.hpp b/src/Decomposition/nn_processor.hpp index a8eff91ef0367e663d7da9189beb5356afdc9552..a3fdb27a43cb0967bbf1abfa210056419e64375c 100644 --- a/src/Decomposition/nn_processor.hpp +++ b/src/Decomposition/nn_processor.hpp @@ -86,8 +86,7 @@ class nn_prcs nnpst.pos.add(c); } - /*! \brief In case of periodic boundary conditions we have to add boxes - * at the borders + /*! \brief In case of periodic boundary conditions we replicate the sub-domains at the border * * \param domain Domain box * \param boundary boundary conditions @@ -109,6 +108,7 @@ class nn_prcs if (check_valid(cmbs[j],bc) == false) continue; + // Calculate the sector box Box<dim,T> bp; Point<dim,T> shift; @@ -119,7 +119,7 @@ class nn_prcs case 1: bp.setLow(k,domain.getHigh(k)+ghost.getLow(k)); bp.setHigh(k,domain.getHigh(k)); - shift.get(k) = -domain.getHigh(k); + shift.get(k) = -domain.getHigh(k)+domain.getLow(k); break; case 0: bp.setLow(k,domain.getLow(k)); @@ -129,7 +129,7 @@ class nn_prcs case -1: bp.setLow(k,domain.getLow(k)); bp.setHigh(k,ghost.getHigh(k)); - shift.get(k) = domain.getHigh(k); + shift.get(k) = domain.getHigh(k)-domain.getLow(k); break; } } @@ -211,6 +211,10 @@ public: */ static bool inline check_valid(comb<dim> cmb,const size_t (& bc)[dim]) { + // the combination 0 is not valid + if (cmb.n_zero() == dim) + return false; + for (size_t i = 0 ; i < dim ; i++) { if (bc[i] == NON_PERIODIC && cmb.getComb()[i] != 0) diff --git a/src/Grid/grid_dist_id_iterator_sub.hpp b/src/Grid/grid_dist_id_iterator_sub.hpp index 6c985fc40d3511135062bedb3addcc1564a527e9..244fa2e589faf36f5b46cf65c43de329d0e6a780 100644 --- a/src/Grid/grid_dist_id_iterator_sub.hpp +++ b/src/Grid/grid_dist_id_iterator_sub.hpp @@ -131,7 +131,11 @@ class grid_dist_iterator_sub grid_dist_iterator_sub(const grid_dist_iterator_sub<dim,device_grid> & tmp) :g_c(tmp.g_c),gList(tmp.gList),gdb_ext(gdb_ext),start(tmp.start),stop(tmp.stop) { - a_it.reinitialize(tmp.a_it); + // get the next grid iterator + if (g_c < gList.size()) + { + a_it.reinitialize(tmp.a_it); + } } /*! \brief Constructor of the distributed grid iterator diff --git a/src/Grid/grid_dist_id_unit_test.hpp b/src/Grid/grid_dist_id_unit_test.hpp index d33d76d758ae670f91b67e79303628b1c8b64b37..7c337b386709f45ffe9e92bdfc1b356e23953a77 100644 --- a/src/Grid/grid_dist_id_unit_test.hpp +++ b/src/Grid/grid_dist_id_unit_test.hpp @@ -3,6 +3,7 @@ #include "grid_dist_id.hpp" #include "data_type/scalar.hpp" +#include "data_type/aggregate.hpp" BOOST_AUTO_TEST_SUITE( grid_dist_id_test ) @@ -587,6 +588,113 @@ void Test3D_gg(const Box<3,float> & domain, long int k, long int gk) } } +/*! \brief Test when the domain is not from 0.0 to 1.0 + * + * + */ + +void Test3D_domain(const Box<3,float> & domain, long int k) +{ + long int big_step = k / 30; + big_step = (big_step == 0)?1:big_step; + long int small_step = 21; + + print_test( "Testing 3D grid k<=",k); + + // 3D test + for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step ) + { + BOOST_TEST_CHECKPOINT( "Testing 3D grid k=" << k ); + + // grid size + size_t sz[3]; + sz[0] = k; + sz[1] = k; + sz[2] = k; + + // factor + float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f); + + // Ghost + Ghost<3,float> g(0.01 / factor); + + // Distributed grid with id decomposition + grid_dist_id<3, float, aggregate<long int,long int>, CartDecomposition<3,float>> g_dist(sz,domain,g); + + // check the consistency of the decomposition + bool val = g_dist.getDecomposition().check_consistency(); + BOOST_REQUIRE_EQUAL(val,true); + + // Grid sm + grid_sm<3,void> info(sz); + + // get the domain iterator + size_t count = 0; + + auto dom = g_dist.getDomainIterator(); + + while (dom.isNext()) + { + auto key = dom.get(); + auto key_g = g_dist.getGKey(key); + + g_dist.template get<0>(key) = count; + g_dist.template get<1>(key) = info.LinId(key_g); + + // Count the point + count++; + + ++dom; + } + + size_t count2 = count; + openfpm::vector<size_t> pnt; + + // Get the total size of the local grids on each processors + // and the total size + Vcluster & v_cl = g_dist.getVC(); + v_cl.sum(count2); + v_cl.allGather(count,pnt); + v_cl.execute(); + size_t 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); + + // Check + BOOST_REQUIRE_EQUAL(count2,(size_t)k*k*k); + + // sync the ghost + g_dist.template ghost_get<0,1>(); + + bool match = true; + + // check that the communication is correctly completed + + auto domg = g_dist.getDomainGhostIterator(); + + // check that the grid with the ghost past store the correct information + while (domg.isNext()) + { + auto key = domg.get(); + auto key_g = g_dist.getGKey(key); + + // In this case the boundary condition are non periodic + if (g_dist.isInside(key_g)) + { + match &= (g_dist.template get<1>(key) == info.LinId(key_g))?true:false; + } + + ++domg; + } + + BOOST_REQUIRE_EQUAL(match,true); + } +} + + + void Test2D_complex(const Box<2,float> & domain, long int k) { typedef Point_test<float> p; @@ -1253,6 +1361,20 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_decomposition_iterator ) Test3D_decit(domain3,k); } + +BOOST_AUTO_TEST_CASE( grid_dist_id_domain_test_use) +{ + // Initialize the global VCluster + init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv); + + // Domain + Box<3,float> domain3({0.1,0.1,0.1},{1.1,1.1,1.1}); + + long int k = 128*128*128*global_v_cluster->getProcessingUnits(); + k = std::pow(k, 1/3.); + Test3D_domain(domain3,k); +} + BOOST_AUTO_TEST_SUITE_END() #endif diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp index 2090ca4c7d830e32be53c6163e4429c2a3d58108..44caa10319393c7531fe0ab7712800703bb96ecb 100644 --- a/src/Vector/vector_dist.hpp +++ b/src/Vector/vector_dist.hpp @@ -1046,7 +1046,7 @@ public: // Calculate the division array and the cell box for (size_t i = 0 ; i < dim ; i++) { - div[i] = (pbox.getP2().get(i) - pbox.getP1().get(i))/ r_cut; + div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i))/ r_cut); div[i]++; cell_box.setLow(i,0.0); diff --git a/src/Vector/vector_dist_key.hpp b/src/Vector/vector_dist_key.hpp index 03a0286ddd9cf0460332e6b87c538f776952cd4a..80f77580f6b0306c392ecde39a385e8b79d6b0a8 100644 --- a/src/Vector/vector_dist_key.hpp +++ b/src/Vector/vector_dist_key.hpp @@ -24,12 +24,22 @@ class vect_dist_key_dx public: + /*! \brief set the key + * + * \return the local key + * + */ + inline void setKey(size_t key) + { + this->key = key; + } + /*! \brief Get the key * * \return the local key * */ - size_t getKey() + inline size_t getKey() { return key; } @@ -46,10 +56,14 @@ public: return ts.str(); } - vect_dist_key_dx(size_t key) + inline vect_dist_key_dx(size_t key) :key(key) { } + + inline vect_dist_key_dx() + { + } };