Merge pull request #644 from Xiangyu-Hu/wip/kong/split_dynamics_restr…

…ucture

Remove split CCL

src/shared/common/sphinxsys_containers.h

@@ -77,8 +77,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

@@ -13,40 +13,18 @@ BaseCellLinkedList::
     BaseCellLinkedList(BaseParticles &base_particles, SPHAdaptation &sph_adaptation)
     : 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,
                                BaseParticles &base_particles, SPHAdaptation &sph_adaptation)
     : BaseCellLinkedList(base_particles, sph_adaptation), Mesh(tentative_bounds, grid_spacing, 2),
-      use_split_cell_lists_(false), cell_offset_list_size_(NumberOfCells() + 1),
+      cell_offset_list_size_(NumberOfCells() + 1),
       index_list_size_(SMAX(base_particles.ParticlesBound(), cell_offset_list_size_)),
       dv_particle_index_(base_particles.registerDiscreteVariableOnly<UnsignedInt>("ParticleIndex", index_list_size_)),
       dv_cell_offset_(base_particles.registerDiscreteVariableOnly<UnsignedInt>("CellOffset", cell_offset_list_size_)),
-      cell_index_lists_(nullptr), cell_data_lists_(nullptr)
+      cell_index_lists_(nullptr), cell_data_lists_(nullptr),
+      number_of_split_cell_lists_(static_cast<size_t>(pow(3, Dimensions)))
 {
     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()
@@ -89,23 +67,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();
@@ -123,11 +84,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

@@ -33,6 +33,7 @@
 #define MESH_CELL_LINKED_LIST_H
 
 #include "base_mesh.h"
+#include "execution_policy.h"
 #include "neighborhood.h"
 
 namespace SPH
@@ -51,20 +52,14 @@ 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(BaseParticles &base_particles, SPHAdaptation &sph_adaptation);
-    virtual ~BaseCellLinkedList() {};
+    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. */
@@ -105,14 +100,6 @@ class NeighborSearch : public Mesh
 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_;
 
     UnsignedInt cell_offset_list_size_;
     UnsignedInt index_list_size_; // at least number_of_cells_pluse_one_
@@ -124,10 +111,11 @@ class CellLinkedList : public BaseCellLinkedList, public Mesh
     ConcurrentIndexVector *cell_index_lists_;
     /** non-concurrent list data rewritten for building neighbor list */
     ListDataVector *cell_data_lists_;
+    /**< number of split cell lists */
+    size_t number_of_split_cell_lists_;
 
     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)
     {
@@ -140,8 +128,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;
@@ -163,6 +149,12 @@ class CellLinkedList : public BaseCellLinkedList, public Mesh
     UnsignedInt getCellOffsetListSize() { return cell_offset_list_size_; };
     DiscreteVariable<UnsignedInt> *getParticleIndex() { return dv_particle_index_; };
     DiscreteVariable<UnsignedInt> *getCellOffset() { return dv_cell_offset_; };
+
+    /** split algorithm */;
+    template <class LocalDynamicsFunction>
+    void particle_for_split(const execution::SequencedPolicy &, const LocalDynamicsFunction &local_dynamics_function);
+    template <class LocalDynamicsFunction>
+    void particle_for_split(const execution::ParallelPolicy &, const LocalDynamicsFunction &local_dynamics_function);
 };
 
 template <>
@@ -172,7 +164,7 @@ class RefinedMesh<CellLinkedList> : public CellLinkedList
     RefinedMesh(BoundingBox tentative_bounds, CellLinkedList &coarse_mesh,
                 BaseParticles &base_particles, SPHAdaptation &sph_adaptation)
         : CellLinkedList(tentative_bounds, 0.5 * coarse_mesh.GridSpacing(),
-                         base_particles, sph_adaptation) {};
+                         base_particles, sph_adaptation){};
 };
 
 /**
@@ -184,16 +176,14 @@ class MultilevelCellLinkedList : public MultilevelMesh<BaseCellLinkedList, CellL
 {
   protected:
     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);
 
   public:
     MultilevelCellLinkedList(BoundingBox tentative_bounds,
                              Real reference_grid_spacing, size_t total_levels,
                              BaseParticles &base_particles, SPHAdaptation &sph_adaptation);
-    virtual ~MultilevelCellLinkedList() {};
+    virtual ~MultilevelCellLinkedList(){};
 
     virtual void UpdateCellLists(BaseParticles &base_particles) override;
     void insertParticleIndex(size_t particle_index, const Vecd &particle_position) override;

src/shared/meshes/cell_linked_list.hpp

@@ -83,4 +83,108 @@ void CellLinkedList::searchNeighborsByParticles(
                  });
 }
 //=================================================================================================//
+template <class LocalDynamicsFunction>
+void CellLinkedList::particle_for_split(const execution::SequencedPolicy &, const LocalDynamicsFunction &local_dynamics_function)
+{
+    // foward sweeping
+    for (size_t k = 0; k < number_of_split_cell_lists_; k++)
+    {
+        // get the corresponding 2D/3D split cell index (m, n)
+        // e.g., for k = 0, split_cell_index = (0,0), for k = 3, split_cell_index = (1,0), etc.
+        const Arrayi split_cell_index = transfer1DtoMeshIndex(3 * Arrayi::Ones(), k);
+        // get the number of cells belonging to the split cell k
+        // i_max = (M - m - 1) / 3 + 1, j_max = (N - n - 1) / 3 + 1
+        // e.g. all_cells = (M,N) = (6, 9), (m, n) = (1, 1), then i_max = 2, j_max = 3
+        const Arrayi all_cells_k = (all_cells_ - split_cell_index - Arrayi::Ones()) / 3 + Arrayi::Ones();
+        const size_t number_of_cells = all_cells_k.prod(); // i_max * j_max
+
+        // looping over all cells in the split cell k
+        for (size_t l = 0; l < number_of_cells; l++)
+        {
+            // get the 2D/3D cell index of the l-th cell in the split cell k
+            // (i , j) = (m + 3 * (l / j_max), n + 3 * l % i_max)
+            // e.g. all_cells = (M,N) = (6, 9), (m, n) = (1, 1), l = 0, then (i, j) = (1, 1)
+            // l = 1, then (i, j) = (1, 4), l = 3, then (i, j) = (4, 1), etc.
+            const Arrayi cell_index = split_cell_index + 3 * transfer1DtoMeshIndex(all_cells_k, l);
+            // get the list of particles in the cell (i, j)
+            const ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
+            // looping over all particles in the cell (i, j)
+            for (const size_t index_i : cell_list)
+            {
+                local_dynamics_function(index_i);
+            }
+        }
+    }
+
+    // backward sweeping
+    for (size_t k = number_of_split_cell_lists_; k != 0; --k)
+    {
+        const Arrayi split_cell_index = transfer1DtoMeshIndex(3 * Arrayi::Ones(), k - 1);
+        const Arrayi all_cells_k = (all_cells_ - split_cell_index - Arrayi::Ones()) / 3 + Arrayi::Ones();
+        const size_t number_of_cells = all_cells_k.prod();
+
+        for (size_t l = 0; l < number_of_cells; l++)
+        {
+            const Arrayi cell_index = split_cell_index + 3 * transfer1DtoMeshIndex(all_cells_k, l);
+            const ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
+            for (size_t i = cell_list.size(); i != 0; --i)
+            {
+                local_dynamics_function(cell_list[i - 1]);
+            }
+        }
+    }
+}
+//=================================================================================================//
+template <class LocalDynamicsFunction>
+void CellLinkedList::particle_for_split(const execution::ParallelPolicy &, const LocalDynamicsFunction &local_dynamics_function)
+{
+    // foward sweeping
+    for (size_t k = 0; k < number_of_split_cell_lists_; k++)
+    {
+        const Arrayi split_cell_index = transfer1DtoMeshIndex(3 * Arrayi::Ones(), k);
+        const Arrayi all_cells_k = (all_cells_ - split_cell_index - Arrayi::Ones()) / 3 + Arrayi::Ones();
+        const size_t number_of_cells = all_cells_k.prod();
+
+        parallel_for(
+            IndexRange(0, number_of_cells),
+            [&](const IndexRange &r)
+            {
+                for (size_t l = r.begin(); l < r.end(); ++l)
+                {
+                    const Arrayi cell_index = split_cell_index + 3 * transfer1DtoMeshIndex(all_cells_k, l);
+                    const ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
+                    for (const size_t index_i : cell_list)
+                    {
+                        local_dynamics_function(index_i);
+                    }
+                }
+            },
+            ap);
+    }
+
+    // backward sweeping
+    for (size_t k = number_of_split_cell_lists_; k != 0; --k)
+    {
+        const Arrayi split_cell_index = transfer1DtoMeshIndex(3 * Arrayi::Ones(), k - 1);
+        const Arrayi all_cells_k = (all_cells_ - split_cell_index - Arrayi::Ones()) / 3 + Arrayi::Ones();
+        const size_t number_of_cells = all_cells_k.prod();
+
+        parallel_for(
+            IndexRange(0, number_of_cells),
+            [&](const IndexRange &r)
+            {
+                for (size_t l = r.begin(); l < r.end(); ++l)
+                {
+                    const Arrayi cell_index = split_cell_index + 3 * transfer1DtoMeshIndex(all_cells_k, l);
+                    const ConcurrentIndexVector &cell_list = getCellDataList(cell_index_lists_, cell_index);
+                    for (size_t i = cell_list.size(); i != 0; --i)
+                    {
+                        local_dynamics_function(cell_list[i - 1]);
+                    }
+                }
+            },
+            ap);
+    }
+}
+//=================================================================================================//
 } // namespace SPH

src/shared/particle_dynamics/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,28 +199,25 @@ 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)
+    explicit 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");
     };
-    virtual ~InteractionSplit(){};
 
     /** run the main interaction step between particles. */
-    virtual void runMainStep(Real dt) override
+    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(ExecutionPolicy(), [&](size_t i)
+                                             { this->interaction(i, dt * 0.5); });
     }
 };