Remove split CCL

src/for_2D_build/meshes/mesh_iterators.hpp

@@ -84,4 +84,46 @@ void mesh_parallel_for(const MeshRange &mesh_range, const LocalFunction &local_f
         ap);
 }
 //=================================================================================================//
+template <typename LocalFunction, typename... Args>
+void mesh_split_parallel_for(const MeshRange &mesh_range, const Arrayi &stride, const LocalFunction &local_function, Args &&...args)
+{
+    // forward sweeping
+    for (int m = 0; m < stride[0]; m++)
+        for (int n = 0; n < stride[1]; n++)
+        {
+            parallel_for(
+                IndexRange2d((mesh_range.first)[0] + m, (mesh_range.second)[0],
+                             (mesh_range.first)[1] + n, (mesh_range.second)[1]),
+                [&](const IndexRange2d &r)
+                {
+                    for (size_t i = r.rows().begin(); i != r.rows().end(); i++)
+                        for (size_t j = r.cols().begin(); j != r.cols().end(); j++)
+                        {
+                            if ((i - m) % stride[0] == 0 && (j - n) % stride[1] == 0)
+                                local_function(Array2i(i, j));
+                        }
+                },
+                ap);
+        }
+
+    // backward sweeping
+    for (int m = stride[0] - 1; m >= 0; m--)
+        for (int n = stride[1] - 1; n >= 0; n--)
+        {
+            parallel_for(
+                IndexRange2d((mesh_range.first)[0] + m, (mesh_range.second)[0],
+                             (mesh_range.first)[1] + n, (mesh_range.second)[1]),
+                [&](const IndexRange2d &r)
+                {
+                    for (size_t i = r.rows().end(); i != r.rows().begin(); i--)
+                        for (size_t j = r.cols().end(); j != r.cols().begin(); j--)
+                        {
+                            if ((i - 1 - m) % stride[0] == 0 && (j - 1 - n) % stride[1] == 0)
+                                local_function(Array2i(i - 1, j - 1));
+                        }
+                },
+                ap);
+        }
+}
+//=================================================================================================//
 } // namespace SPH

src/for_3D_build/meshes/mesh_iterators.hpp

@@ -94,4 +94,52 @@ void mesh_parallel_for(const MeshRange &mesh_range, const LocalFunction &local_f
         ap);
 }
 //=================================================================================================//
+template <typename LocalFunction, typename... Args>
+void mesh_split_parallel_for(const MeshRange &mesh_range, const Arrayi &stride, const LocalFunction &local_function, Args &&...args)
+{
+    // forward sweeping
+    for (int m = 0; m < stride[0]; m++)
+        for (int n = 0; n < stride[1]; n++)
+            for (int p = 0; p < stride[2]; p++)
+            {
+                parallel_for(
+                    IndexRange3d((mesh_range.first)[0] + m, (mesh_range.second)[0],
+                                 (mesh_range.first)[1] + n, (mesh_range.second)[1],
+                                 (mesh_range.first)[2] + p, (mesh_range.second)[2]),
+                    [&](const IndexRange3d &r)
+                    {
+                        for (size_t i = r.pages().begin(); i != r.pages().end(); ++i)
+                            for (size_t j = r.rows().begin(); j != r.rows().end(); ++j)
+                                for (size_t k = r.cols().begin(); k != r.cols().end(); ++k)
+                                {
+                                    if ((i - m) % stride[0] == 0 && (j - n) % stride[1] == 0 && (k - p) % stride[2] == 0)
+                                        local_function(Array3i(i, j, k));
+                                }
+                    },
+                    ap);
+            }
+
+    // backward sweeping
+    for (int m = stride[0] - 1; m >= 0; m--)
+        for (int n = stride[1] - 1; n >= 0; n--)
+            for (int p = stride[2] - 1; p >= 0; p--)
+            {
+                parallel_for(
+                    IndexRange3d((mesh_range.first)[0] + m, (mesh_range.second)[0],
+                                 (mesh_range.first)[1] + n, (mesh_range.second)[1],
+                                 (mesh_range.first)[2] + p, (mesh_range.second)[2]),
+                    [&](const IndexRange3d &r)
+                    {
+                        for (size_t i = r.pages().end(); i != r.pages().begin(); i--)
+                            for (size_t j = r.rows().end(); j != r.rows().begin(); j--)
+                                for (size_t k = r.cols().end(); k != r.cols().begin(); k--)
+                                {
+                                    if ((i - 1 - m) % stride[0] == 0 && (j - 1 - n) % stride[1] == 0 && (k - 1 - p) % stride[2] == 0)
+                                        local_function(Array3i(i - 1, j - 1, k - 1));
+                                }
+                    },
+                    ap);
+            }
+}
+//=================================================================================================//
 } // namespace SPH

src/shared/common/sph_data_containers.h

@@ -62,8 +62,6 @@ using ListData = std::pair<size_t, Vecd>;
 using ListDataVector = StdLargeVec<ListData>;
 using DataListsInCells = StdLargeVec<ListDataVector *>;
 using ConcurrentCellLists = ConcurrentVec<ConcurrentIndexVector *>;
-/** Cell list for splitting algorithms. */
-using SplitCellLists = StdVec<ConcurrentCellLists>;
 /** Cell list for periodic boundary condition algorithms. */
 using CellLists = std::pair<ConcurrentCellLists, DataListsInCells>;
 

src/shared/meshes/cell_linked_list.cpp

@@ -14,36 +14,13 @@ BaseCellLinkedList::
     : BaseMeshField("CellLinkedList"),
       kernel_(*sph_adaptation.getKernel()) {}
 //=================================================================================================//
-SplitCellLists *BaseCellLinkedList::getSplitCellLists()
-{
-    std::cout << "\n Error: SplitCellList not defined!" << std::endl;
-    std::cout << __FILE__ << ':' << __LINE__ << std::endl;
-    exit(1);
-    return nullptr;
-}
-//=================================================================================================//
-void BaseCellLinkedList::setUseSplitCellLists()
-{
-    std::cout << "\n Error: SplitCellList not defined!" << std::endl;
-    std::cout << __FILE__ << ':' << __LINE__ << std::endl;
-    exit(1);
-};
-//=================================================================================================//
-void BaseCellLinkedList::clearSplitCellLists(SplitCellLists &split_cell_lists)
-{
-    for (size_t i = 0; i < split_cell_lists.size(); i++)
-        split_cell_lists[i].clear();
-}
-//=================================================================================================//
 CellLinkedList::CellLinkedList(BoundingBox tentative_bounds, Real grid_spacing,
                                SPHAdaptation &sph_adaptation)
     : BaseCellLinkedList(sph_adaptation), Mesh(tentative_bounds, grid_spacing, 2),
-      use_split_cell_lists_(false), cell_index_lists_(nullptr), cell_data_lists_(nullptr)
+      cell_index_lists_(nullptr), cell_data_lists_(nullptr)
 {
     allocateMeshDataMatrix();
     single_cell_linked_list_level_.push_back(this);
-    size_t number_of_split_cell_lists = pow(3, Dimensions);
-    split_cell_lists_.resize(number_of_split_cell_lists);
 }
 //=================================================================================================//
 void CellLinkedList ::allocateMeshDataMatrix()
@@ -86,23 +63,6 @@ void CellLinkedList::UpdateCellListData(BaseParticles &base_particles)
         });
 }
 //=================================================================================================//
-void CellLinkedList::updateSplitCellLists(SplitCellLists &split_cell_lists)
-{
-    clearSplitCellLists(split_cell_lists);
-    mesh_parallel_for(
-        MeshRange(Arrayi::Zero(), all_cells_),
-        [&](const Arrayi &cell_index)
-        {
-            ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
-            size_t real_particles_in_cell = cell_list.size();
-            if (real_particles_in_cell != 0)
-            {
-                split_cell_lists[transferMeshIndexTo1D(3 * Arrayi::Ones(), mod(cell_index, 3))]
-                    .push_back(&cell_list);
-            }
-        });
-}
-//=================================================================================================//
 void CellLinkedList::UpdateCellLists(BaseParticles &base_particles)
 {
     clearCellLists();
@@ -120,11 +80,6 @@ void CellLinkedList::UpdateCellLists(BaseParticles &base_particles)
         ap);
 
     UpdateCellListData(base_particles);
-
-    if (use_split_cell_lists_)
-    {
-        updateSplitCellLists(split_cell_lists_);
-    }
 }
 //=================================================================================================//
 void CellLinkedList ::insertParticleIndex(size_t particle_index, const Vecd &particle_position)

src/shared/meshes/cell_linked_list.h

@@ -52,20 +52,13 @@ class BaseCellLinkedList : public BaseMeshField
   protected:
     Kernel &kernel_;
 
-    /** clear split cell lists in this mesh*/
-    virtual void clearSplitCellLists(SplitCellLists &split_cell_lists);
-    /** update split particle list in this mesh */
-    virtual void updateSplitCellLists(SplitCellLists &split_cell_lists) = 0;
-
   public:
     BaseCellLinkedList(SPHAdaptation &sph_adaptation);
     virtual ~BaseCellLinkedList(){};
 
     /** access concrete cell linked list levels*/
     virtual StdVec<CellLinkedList *> CellLinkedListLevels() = 0;
     virtual void UpdateCellLists(BaseParticles &base_particles) = 0;
-    virtual SplitCellLists *getSplitCellLists();
-    virtual void setUseSplitCellLists();
     /** Insert a cell-linked_list entry to the concurrent index list. */
     virtual void insertParticleIndex(size_t particle_index, const Vecd &particle_position) = 0;
     /** Insert a cell-linked_list entry of the index and particle position pair. */
@@ -88,14 +81,6 @@ class BaseCellLinkedList : public BaseMeshField
 class CellLinkedList : public BaseCellLinkedList, public Mesh
 {
     StdVec<CellLinkedList *> single_cell_linked_list_level_;
-    /**
-     * @brief particle by cells lists is for parallel splitting algorithm.
-     * All particles in each cell are collected together.
-     * If two particles each belongs two different cell entries,
-     * they have no interaction because they are too far.
-     */
-    SplitCellLists split_cell_lists_;
-    bool use_split_cell_lists_;
 
   protected:
     /** using concurrent vectors due to writing conflicts when building the list */
@@ -105,7 +90,6 @@ class CellLinkedList : public BaseCellLinkedList, public Mesh
 
     void allocateMeshDataMatrix(); /**< allocate memories for addresses of data packages. */
     void deleteMeshDataMatrix();   /**< delete memories for addresses of data packages. */
-    virtual void updateSplitCellLists(SplitCellLists &split_cell_lists) override;
     template <typename DataListsType>
     DataListsType &getCellDataList(DataListsType *data_lists, const Arrayi &cell_index)
     {
@@ -117,8 +101,6 @@ class CellLinkedList : public BaseCellLinkedList, public Mesh
     ~CellLinkedList() { deleteMeshDataMatrix(); };
 
     void clearCellLists();
-    virtual SplitCellLists *getSplitCellLists() override { return &split_cell_lists_; };
-    virtual void setUseSplitCellLists() override { use_split_cell_lists_ = true; };
     void UpdateCellListData(BaseParticles &base_particles);
     virtual void UpdateCellLists(BaseParticles &base_particles) override;
     void insertParticleIndex(size_t particle_index, const Vecd &particle_position) override;
@@ -134,6 +116,10 @@ class CellLinkedList : public BaseCellLinkedList, public Mesh
     template <class DynamicsRange, typename GetSearchDepth, typename GetNeighborRelation>
     void searchNeighborsByParticles(DynamicsRange &dynamics_range, ParticleConfiguration &particle_configuration,
                                     GetSearchDepth &get_search_depth, GetNeighborRelation &get_neighbor_relation);
+
+    /** split algorithm */;
+    template <class LocalDynamicsFunction>
+    void particle_for_split(const LocalDynamicsFunction &local_dynamics_function);
 };
 
 template <>
@@ -153,8 +139,6 @@ class MultilevelCellLinkedList : public MultilevelMesh<BaseCellLinkedList, CellL
 {
   protected:
     StdLargeVec<Real> &h_ratio_; /**< Smoothing length for each level. */
-    /** Update split cell list. */
-    virtual void updateSplitCellLists(SplitCellLists &split_cell_lists) override{};
     /** determine mesh level from particle cutoff radius */
     inline size_t getMeshLevel(Real particle_cutoff_radius);
 

src/shared/meshes/cell_linked_list.hpp

@@ -44,4 +44,20 @@ void CellLinkedList::searchNeighborsByParticles(
                  });
 }
 //=================================================================================================//
+template <class LocalDynamicsFunction>
+void CellLinkedList::particle_for_split(const LocalDynamicsFunction &local_dynamics_function)
+{
+    mesh_split_parallel_for(
+        MeshRange(Arrayi::Zero(), all_cells_),
+        3 * Arrayi::Ones(),
+        [&](const Arrayi &cell_index)
+        {
+            const ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
+            for (size_t index_i : cell_list)
+            {
+                local_dynamics_function(index_i);
+            }
+        });
+}
+//=================================================================================================//
 } // namespace SPH

src/shared/meshes/mesh_iterators.h

@@ -107,5 +107,8 @@ void mesh_for(const MeshRange &mesh_range, const LocalFunction &local_function,
 /** Iterator on the mesh by looping index. parallel computing. */
 template <typename LocalFunction, typename... Args>
 void mesh_parallel_for(const MeshRange &mesh_range, const LocalFunction &local_function, Args &&...args);
+/** Iterator on the mesh by looping index with a stride. parallel computing. */
+template <typename LocalFunction, typename... Args>
+void mesh_split_parallel_for(const MeshRange &mesh_range, const Arrayi &stride, const LocalFunction &local_function, Args &&...args);
 } // namespace SPH
 #endif // MESH_ITERATORS_H

src/shared/particle_dynamics/particle_dynamics_algorithms.h

@@ -52,6 +52,7 @@
 
 #include "base_local_dynamics.h"
 #include "base_particle_dynamics.h"
+#include "cell_linked_list.hpp"
 #include "particle_iterators.h"
 
 #include <type_traits>
@@ -198,16 +199,15 @@ class InteractionSplit : public BaseInteractionDynamics<LocalDynamicsType, Paral
 {
   protected:
     RealBody &real_body_;
-    SplitCellLists &split_cell_lists_;
+    CellLinkedList &cell_linked_list_;
 
   public:
     template <typename... Args>
     InteractionSplit(Args &&... args)
         : BaseInteractionDynamics<LocalDynamicsType, ParallelPolicy>(std::forward<Args>(args)...),
           real_body_(DynamicCast<RealBody>(this, this->getSPHBody())),
-          split_cell_lists_(*real_body_.getCellLinkedList().getSplitCellLists())
+          cell_linked_list_(DynamicCast<CellLinkedList>(this, real_body_.getCellLinkedList()))
     {
-        real_body_.getCellLinkedList().setUseSplitCellLists();
         static_assert(!has_initialize<LocalDynamicsType>::value &&
                           !has_update<LocalDynamicsType>::value,
                       "LocalDynamicsType does not fulfill InteractionSplit requirements");
@@ -217,9 +217,8 @@ class InteractionSplit : public BaseInteractionDynamics<LocalDynamicsType, Paral
     /** run the main interaction step between particles. */
     virtual void runMainStep(Real dt) override
     {
-        particle_for(ExecutionPolicy(),
-                     split_cell_lists_,
-                     [&](size_t i) { this->interaction(i, dt * 0.5); });
+        cell_linked_list_.particle_for_split([&](size_t i)
+                                             { this->interaction(i, dt * 0.5); });
     }
 };