[WIP] Fractal with bconds

libmpdata++-config.cmake

@@ -32,7 +32,7 @@ set(libmpdataxx_INCLUDE_DIRS "${CMAKE_CURRENT_LIST_DIR}/../../include/")
 
 ############################################################################################
 # debug mode compiler flags
-set(libmpdataxx_CXX_FLAGS_DEBUG "${libmpdataxx_CXX_FLAGS_DEBUG} -std=c++14 -DBZ_DEBUG -g -Wno-enum-compare -Wfatal-errors") #TODO: -Og if compiler supports it?
+set(libmpdataxx_CXX_FLAGS_DEBUG "${libmpdataxx_CXX_FLAGS_DEBUG} -std=c++17 -DBZ_DEBUG -g -Wno-enum-compare -Wfatal-errors") #TODO: -Og if compiler supports it?
 
 
 ############################################################################################
@@ -42,7 +42,7 @@ if(
   CMAKE_CXX_COMPILER_ID STREQUAL "Clang" OR
   CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang" 
 )
-  set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=c++14 -DNDEBUG -Ofast -march=native -Wno-enum-compare")
+  set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=c++17 -DNDEBUG -Ofast -march=native -Wno-enum-compare")
 
   # preventing Kahan summation from being optimised out
   if (
@@ -58,25 +58,25 @@ if(
   CMAKE_CXX_COMPILER_ID STREQUAL "Intel"
 )
   # flags taken from -fast but without -static
-  set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=gnu++14 -DNDEBUG -xHOST -O3 -ipo -no-prec-div -fp-model fast=2")
+  set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=gnu++17 -DNDEBUG -xHOST -O3 -ipo -no-prec-div -fp-model fast=2")
 endif()
 
 
 ############################################################################################
 # C++14
-include(CheckCXXSourceCompiles)
-set(CMAKE_REQUIRED_FLAGS "-std=c++14")
-check_cxx_source_compiles("
-  #include <type_traits>
-  auto f() { return 1;}
-  template <bool a, class b> using ei=std::enable_if<a,b>; 
-  struct a {a(int){}};struct b:a {using a::a;};  
-  int main(){b i(1);}
-" CXX14_SUPPORTED)
-if (NOT CXX14_SUPPORTED)
-  message(FATAL_ERROR "C++14 compatibility test failed - please update your compiler or point CMake to another one with -DCMAKE_CXX_COMPILER=...")
-endif()
-unset(CMAKE_REQUIRED_FLAGS)
+#include(CheckCXXSourceCompiles)
+#set(CMAKE_REQUIRED_FLAGS "-std=c++14")
+#check_cxx_source_compiles("
+#  #include <type_traits>
+#  auto f() { return 1;}
+#  template <bool a, class b> using ei=std::enable_if<a,b>; 
+#  struct a {a(int){}};struct b:a {using a::a;};  
+#  int main(){b i(1);}
+#" CXX14_SUPPORTED)
+#if (NOT CXX14_SUPPORTED)
+#  message(FATAL_ERROR "C++14 compatibility test failed - please update your compiler or point CMake to another one with -DCMAKE_CXX_COMPILER=...")
+#endif()
+#unset(CMAKE_REQUIRED_FLAGS)
 
 
 ############################################################################################

libmpdata++/bcond/detail/remote_3d_common.hpp

@@ -32,10 +32,8 @@ namespace libmpdataxx
         // based on the difference between idx to be sent by this thread and idx of this process
         idx_t extend_idx(idx_t idx)
         {
-//std::cerr << "extend idx start idx(1): " << idx.lbound(1) << ", " << idx.ubound(1) << std::endl;
           idx.lbound(1) = 0                 + idx.lbound(1) - thread_j.first(); // does it have to start at 0?
           idx.ubound(1) = (grid_size_y - 1) + idx.ubound(1) - thread_j.last();  // does it have to end at grid_size_y - 1?
-//std::cerr << "extend idx end idx(1): " << idx.lbound(1) << ", " << idx.ubound(1) << std::endl;
           return idx;
         }
 
@@ -68,13 +66,14 @@ namespace libmpdataxx
 
         // ctor
         remote_3d_common(
+          boost::mpi::communicator &mpic,
           const rng_t &i,
           const std::array<int, 3> &distmem_grid_size,
           const rng_t _thread_j,
           const int thread_rank,
           const int thread_size
         ) :
-          parent_t(i, distmem_grid_size, true, thread_rank, thread_size), // true indicating that this is a bcond done with a single thread
+          parent_t(mpic, i, distmem_grid_size, true, thread_rank, thread_size), // true indicating that this is a bcond done with a single thread
           thread_j(_thread_j),
           grid_size_y(distmem_grid_size[1])
         {

libmpdata++/bcond/detail/remote_common.hpp

@@ -78,12 +78,6 @@ namespace libmpdataxx
           const int
             msg_send = dir == left ? left : rght;
 
-//          std::cerr << "send_hlpr idx dir " << dir << " : " 
-//            << " (" << idx_send.lbound(0) << ", " << idx_send.ubound(0) << ")"  
-//            << " (" << idx_send.lbound(1) << ", " << idx_send.ubound(1) << ")"  
-//            << " (" << idx_send.lbound(2) << ", " << idx_send.ubound(2) << ")"  
-//            << std::endl;
-
           // arr_send references part of the send buffer that will be used
           arr_t arr_send(buf_send, a(idx_send).shape(), blitz::neverDeleteData);
           // copying data to be sent
@@ -118,13 +112,6 @@ namespace libmpdataxx
           const int
             msg_recv = dir == left ? rght : left;
 
-//          std::cerr << "recv_hlpr idx dir " << dir << " : " 
-//            << " (" << idx_recv.lbound(0) << ", " << idx_recv.ubound(0) << ")"  
-//            << " (" << idx_recv.lbound(1) << ", " << idx_recv.ubound(1) << ")"  
-//            << " (" << idx_recv.lbound(2) << ", " << idx_recv.ubound(2) << ")"  
-//            << std::endl;
-
-
           // launching async data transfer
           if(a(idx_recv).size()!=0) // TODO: test directly size of idx_recv
           {
@@ -218,24 +205,20 @@ namespace libmpdataxx
 
         // ctor
         remote_common(
+          boost::mpi::communicator &mpic,
           const rng_t &i,
           const std::array<int, n_dims> &distmem_grid_size,
           bool single_threaded = false,
           const int thread_rank = -1, 
-          const int thread_size = -1 
+          const int thread_size = -1
         ) :
           parent_t(i, distmem_grid_size, single_threaded, thread_rank, thread_size)
+          //,mpicom(mpic, boost::mpi::comm_take_ownership)
+          ,mpicom(mpic, boost::mpi::comm_attach)
         {
 #if defined(USE_MPI)
 
           const int slice_size = n_dims==1 ? 1 : (n_dims==2? distmem_grid_size[1]+6 : (distmem_grid_size[1]+6) * (distmem_grid_size[2]+6) ); // 3 is the max halo size (?), so 6 on both sides
-//std::cerr << "remote_common ctor, " 
-//  << " distmem_grid_size[0]: " << distmem_grid_size[0]
-//  << " distmem_grid_size[1]: " << distmem_grid_size[1]
-//  << " distmem_grid_size[2]: " << distmem_grid_size[2]
-//  << " slice_size: " << slice_size 
-//  << " halo: " << halo
-//  << std::endl;
           // allocate enough memory in buffers to store largest halos to be sent
           buf_send = (real_t *) malloc(halo * slice_size * sizeof(real_t));
           buf_recv = (real_t *) malloc(halo * slice_size * sizeof(real_t));

libmpdata++/blitz.hpp

@@ -51,6 +51,7 @@
 namespace libmpdataxx
 {
   template <int n_dims> using idx_t = blitz::RectDomain<n_dims>;
+  template <int n_dims> using idxs_t = blitz::StridedDomain<n_dims>;
   using rng_t = blitz::Range;
 
   // non-int ix_t means either rng_t or idx_t

libmpdata++/concurr/boost_thread.hpp

@@ -51,10 +51,8 @@ namespace libmpdataxx
 
 
         // ctor
-        mem_t(const std::array<int, solver_t::n_dims> &grid_size) :
-          b(size(grid_size[0])),
-          parent_t::mem_t(grid_size, size(grid_size[0]))
-        {};
+        mem_t(const std::array<int, solver_t::n_dims> &grid_size, const int n_ref) : b(size(grid_size[0])), parent_t::mem_t(grid_size, size(grid_size[0]), n_ref) {};
+        mem_t(const std::array<int, solver_t::n_dims> &grid_size)                  : b(size(grid_size[0])), parent_t::mem_t(grid_size, size(grid_size[0]))        {};
 
         void barrier()
         {
@@ -81,7 +79,7 @@ namespace libmpdataxx
 
       // ctor
       boost_thread(const typename solver_t::rt_params_t &p) :
-        parent_t(p, new mem_t(p.grid_size), mem_t::size(p.grid_size[solver_t::n_dims < 3 ? 0 : 1])) // note 3D domain decomposition in y direction
+        parent_t(p, detail::mem_factory<mem_t, solver_t>(p), mem_t::size(p.grid_size[solver_t::n_dims < 3 ? 0 : 1])) // note 3D domain decomposition in y direction
       {}
 
     };

libmpdata++/concurr/cxx11_thread.hpp

@@ -92,10 +92,8 @@ namespace libmpdataxx
         }
 
         // ctor
-        mem_t(const std::array<int, solver_t::n_dims> &grid_size) :
-          b(size(grid_size[0])),
-          parent_t::mem_t(grid_size, size(grid_size[0]))
-        {};
+        mem_t(const std::array<int, solver_t::n_dims> &grid_size, const int n_ref) : b(size(grid_size[0])), parent_t::mem_t(grid_size, size(grid_size[0]), n_ref) {};
+        mem_t(const std::array<int, solver_t::n_dims> &grid_size)                  : b(size(grid_size[0])), parent_t::mem_t(grid_size, size(grid_size[0]))        {};
 
         void barrier()
         {
@@ -119,7 +117,7 @@ namespace libmpdataxx
 
       // ctor
       cxx11_thread(const typename solver_t::rt_params_t &p) :
-        parent_t(p, new mem_t(p.grid_size), mem_t::size(p.grid_size[solver_t::n_dims < 3 ? 0 : 1])) // note 3D domain decomposition in y direction
+        parent_t(p, detail::mem_factory<mem_t, solver_t>(p), mem_t::size(p.grid_size[solver_t::n_dims < 3 ? 0 : 1])) // note 3D domain decomposition in y direction
       {}
 
     };