Merge pull request #896 from Xiangyu-Hu/xiangyu/unifed_kernel

Xiangyu/unified kernel

Author: Xiangyu-Hu

README.md

@@ -1,28 +1,24 @@
 # ![SPHinXsys Logo](assets/logo.png) SPHinXsys
 
-## Heterogeneous parallelism for mesh dynamics and particle relaxation
-
-It is known that SPHinXsys relies on level-set technique to realize
-particle generation and relaxation for later SPH simulations.
-The level-set field is build on a multi-level Cartesian mesh
-with sparse storage and adaptive mesh refinement (AMR).
-After the level-set field is constructed,
-the body-fitted initial particle distribution,
-especially for solid bodies,
-is obtained by a physical-driving relaxation process
-using level-set field for body-fitting.  
-
-Now, the level-set based operations or more generally mesh dynamics,
-such as small-feature cleaning, re-initialization
-and normal direction computing,
-and the particle relaxation process now can be carried on either CPU or GPU
-with the same numerical algorithm and the help of SYCL kernels.
-
-The test utilizes heterogeneous parallelism for mesh dynamics have been added to the cases
- `test_2d_particle_generator_single_resolution_sycl`,`test_3d_particle_relaxation_single_resolution_sycl`
-and `taylor_bar_sycl` in `test/test_sycl` folder,
-
-showcasing the capabilities of our specially designed framework.
+## SPHinXsys new feature: turning “mesh nightmares” into smooth simulations
+
+Ever had a perfect-looking CAD model that refused to mesh?
+
+Hidden leaks (those tiny yellow marks in the figure) can completely derail numerical simulations — and repairing them can be tedious and frustrating.
+
+![Leaking Geometry](https://media.licdn.com/dms/image/v2/D4D22AQF5CCANH4NQgQ/feedshare-shrink_2048_1536/B4DZiQ_oVfHYAw-/0/1754779257657?e=1758153600&v=beta&t=ugXbUmTgQQ4s7bjkO_f-zY7Y5yIOsUzlE1-sobnqLWc)
+
+With SPHinXsys level set techniques, you can skip the repair work.  
+We directly convert your original (even leaky) geometry into particles, then apply a physics-based relaxation to produce a clean, body-fitted, and uniformly distributed particle set — ready for simulation.
+
+![Levelset and Particles](https://media.licdn.com/dms/image/v2/D4D22AQFmcHYYO-ji2Q/feedshare-shrink_2048_1536/B4DZiQ_oWQGsAo-/0/1754779257939?e=1758153600&v=beta&t=_IWXfhpGhLDzw4ukFiOrrOJlOx3jtS6J48xk5J9pPC4)
+
+No more chasing tiny gaps. Just upload your model and let the algorithm do the work.
+
+You can check the test cases
+`test_3d_particle_relaxation_single_resolution_sycl`
+for different geometries in the `data` folder,
+which showcase the capabilities of our specially designed framework.
 
 Your tests, comments and modification of these test cases would be very welcomed!
 

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-file(GLOB_RECURSE SPHINXSYS_SHARED_HEADERS CONFIGURE_DEPENDS shared/*.h shared/*.hpp)
+file(GLOB_RECURSE SPHINXSYS_SHARED_HEADERS CONFIGURE_DEPENDS shared/*.h shared/*.hpp  shared/*.hxx)
 file(GLOB_RECURSE SPHINXSYS_SHARED_SOURCES CONFIGURE_DEPENDS shared/*.cpp)
 
 if(SPHINXSYS_USE_SYCL)

src/for_2D_build/mesh_dynamics/mesh_local_dynamics.hpp

@@ -134,7 +134,7 @@ inline Real UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->W(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.W(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_);
                 }
             });
@@ -163,7 +163,7 @@ inline Vecd UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->dW(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.dW(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_) *
                                     displacement / (distance + TinyReal);
                 }
@@ -194,7 +194,7 @@ inline Matd UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->d2W(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.d2W(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_) *
                                     displacement * displacement.transpose() / (distance * distance + TinyReal);
                 }

src/for_2D_build/meshes/sparse_storage_mesh/mesh_with_data_packages.hpp

@@ -1,7 +1,7 @@
 #ifndef MESH_WITH_DATA_PACKAGES_2D_HPP
 #define MESH_WITH_DATA_PACKAGES_2D_HPP
 
-#include "mesh_with_data_packages.h"
+#include "mesh_with_data_packages.hxx"
 
 namespace SPH
 {

src/for_3D_build/mesh_dynamics/mesh_local_dynamics.hpp

@@ -146,7 +146,7 @@ inline Real UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j, k).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->W(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.W(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_);
                 }
             });
@@ -181,7 +181,7 @@ inline Vecd UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j, k).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->dW(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.dW(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_) *
                                     displacement / (distance + TinyReal);
                 }
@@ -218,7 +218,7 @@ inline Matd UpdateKernelIntegrals::UpdateKernel::
                     Vecd displacement = -Arrayi(i, j, k).cast<Real>().matrix() * data_spacing_;
                     Real distance = displacement.norm();
                     if (distance < cutoff_radius_)
-                        integral += kernel_->d2W(global_h_ratio_, distance, displacement) *
+                        integral += kernel_.d2W(displacement) *
                                     CutCellVolumeFraction(phi_neighbor, phi_gradient, data_spacing_) *
                                     displacement * displacement.transpose() / (distance * distance + TinyReal);
                 }

src/for_3D_build/meshes/sparse_storage_mesh/mesh_with_data_packages.hpp

@@ -1,7 +1,7 @@
 #ifndef MESH_WITH_DATA_PACKAGES_3D_HPP
 #define MESH_WITH_DATA_PACKAGES_3D_HPP
 
-#include "mesh_with_data_packages.h"
+#include "mesh_with_data_packages.hxx"
 
 #include "block_cluster_3d.h"
 

src/shared/common/base_dynamics.h

@@ -0,0 +1,56 @@
+/* ------------------------------------------------------------------------- *
+ *                                SPHinXsys                                  *
+ * ------------------------------------------------------------------------- *
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for    *
+ * physical accurate simulation and aims to model coupled industrial dynamic *
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH   *
+ * (smoothed particle hydrodynamics), a meshless computational method using  *
+ * particle discretization.                                                  *
+ *                                                                           *
+ * SPHinXsys is partially funded by German Research Foundation               *
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,            *
+ *  HU1527/12-1 and HU1527/12-4.                                             *
+ *                                                                           *
+ * Portions copyright (c) 2017-2025 Technical University of Munich and       *
+ * the authors' affiliations.                                                *
+ *                                                                           *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may   *
+ * not use this file except in compliance with the License. You may obtain a *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+/**
+ * @file 	base_dynamics.h
+ * @brief 	This is for the base classes of all dynamics.
+ * @author	Xiangyu Hu
+ */
+#ifndef BASE_DYNAMICS_H
+#define BASE_DYNAMICS_H
+
+namespace SPH
+{
+template <class ReturnType = void>
+class BaseDynamics
+{
+  public:
+    BaseDynamics() : is_newly_updated_(false) {};
+    virtual ~BaseDynamics() {};
+    bool checkNewlyUpdated() { return is_newly_updated_; };
+    void setNotNewlyUpdated() { is_newly_updated_ = false; };
+
+    template <class IdentifierType>
+    void setUpdated(IdentifierType &identifier)
+    {
+        identifier.setNewlyUpdated();
+        is_newly_updated_ = true;
+    };
+
+    /** There is the interface functions for computing. */
+    virtual ReturnType exec(Real dt = 0.0) = 0;
+
+  private:
+    bool is_newly_updated_;
+};
+} // namespace SPH
+#endif // BASE_DYNAMICS_H

src/shared/geometries/level_set.cpp

@@ -11,43 +11,51 @@ MultilevelLevelSet::MultilevelLevelSet(
     Shape &shape, SPHAdaptation &sph_adaptation, Real refinement_ratio)
     : BaseMeshField("LevelSet_" + shape.getName()), shape_(shape), total_levels_(1)
 {
-    Real reference_data_spacing = coarse_data->DataSpacing() * 0.5;
-    Real global_h_ratio = sph_adaptation.ReferenceSpacing() / reference_data_spacing / refinement_ratio;
-    kernel_ = makeUnique<SingularVariable<KernelTabulatedCK>>(
-        "levelset_kernel", KernelTabulatedCK(*sph_adaptation.getKernel()));
+    Real data_spacing = coarse_data->DataSpacing() * 0.5;
+    Real global_h_ratio = sph_adaptation.ReferenceSpacing() / data_spacing / refinement_ratio;
+    Real smoothing_length = sph_adaptation.ReferenceSmoothingLength() / global_h_ratio;
     global_h_ratio_vec_.push_back(global_h_ratio);
+    neighbor_method_set_.push_back(
+        neighbor_method_keeper_.template createPtr<NeighborMethod<SingleValued>>(
+            *sph_adaptation.getKernel(), smoothing_length, data_spacing));
 
-    initializeLevel(reference_data_spacing, tentative_bounds, coarse_data);
+    initializeLevel(data_spacing, tentative_bounds, coarse_data);
 }
 //=================================================================================================//
 MultilevelLevelSet::MultilevelLevelSet(
-    BoundingBox tentative_bounds, Real reference_data_spacing,
+    BoundingBox tentative_bounds, Real data_spacing,
     size_t total_levels, Shape &shape, SPHAdaptation &sph_adaptation, Real refinement_ratio)
     : BaseMeshField("LevelSet_" + shape.getName()), shape_(shape), total_levels_(total_levels)
 {
-    Real global_h_ratio = sph_adaptation.ReferenceSpacing() / reference_data_spacing / refinement_ratio;
+    Real global_h_ratio = sph_adaptation.ReferenceSpacing() / data_spacing / refinement_ratio;
+    Real smoothing_length = sph_adaptation.ReferenceSmoothingLength() / global_h_ratio;
     global_h_ratio_vec_.push_back(global_h_ratio);
-    kernel_ = makeUnique<SingularVariable<KernelTabulatedCK>>(
-        "levelset_kernel", KernelTabulatedCK(*sph_adaptation.getKernel()));
+    neighbor_method_set_.push_back(
+        neighbor_method_keeper_.template createPtr<NeighborMethod<SingleValued>>(
+            *sph_adaptation.getKernel(), smoothing_length, data_spacing));
 
-    initializeLevel(reference_data_spacing, tentative_bounds);
+    initializeLevel(data_spacing, tentative_bounds);
     for (size_t level = 1; level < total_levels_; ++level)
     {
-        reference_data_spacing *= 0.5; // Halve the data spacing
-        global_h_ratio *= 2;           // Double the ratio
+        data_spacing *= 0.5;     // Halve the data spacing
+        global_h_ratio *= 2;     // Double the ratio
+        smoothing_length *= 0.5; // Halve the smoothing length
         global_h_ratio_vec_.push_back(global_h_ratio);
+        neighbor_method_set_.push_back(
+            neighbor_method_keeper_.template createPtr<NeighborMethod<SingleValued>>(
+                *sph_adaptation.getKernel(), smoothing_length, data_spacing));
 
-        initializeLevel(reference_data_spacing, tentative_bounds, mesh_data_set_[level - 1]);
+        initializeLevel(data_spacing, tentative_bounds, mesh_data_set_[level - 1]);
     }
 }
 //=================================================================================================//
 void MultilevelLevelSet::initializeLevel(
-    Real reference_data_spacing, BoundingBox tentative_bounds, MeshWithGridDataPackagesType *coarse_data)
+    Real data_spacing, BoundingBox tentative_bounds, MeshWithGridDataPackagesType *coarse_data)
 {
     MeshWithGridDataPackagesType *mesh_data =
         mesh_data_ptr_vector_keeper_
             .template createPtr<MeshWithGridDataPackagesType>(
-                tentative_bounds, reference_data_spacing, 4);
+                tentative_bounds, data_spacing, 4);
     mesh_data_set_.push_back(mesh_data);
 
     if (coarse_data == nullptr)

src/shared/geometries/level_set.h

@@ -121,17 +121,18 @@ class MultilevelLevelSet : public BaseMeshField
     void initializeLevel(Real reference_data_spacing, BoundingBox tentative_bounds,
                          MeshWithGridDataPackagesType *coarse_data = nullptr);
     template <class ExecutionPolicy>
-    void initializeMeshVariables(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel);
+    void initializeMeshVariables(const ExecutionPolicy &ex_policy);
     template <class ExecutionPolicy>
-    void initializeKernelIntegralVariables(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel);
+    void initializeKernelIntegralVariables(const ExecutionPolicy &ex_policy);
     template <class ExecutionPolicy>
     void registerProbes(const ExecutionPolicy &ex_policy, size_t level);
 
     Shape &shape_;        /**< the geometry is described by the level set. */
     size_t total_levels_; /**< level 0 is the coarsest */
     StdVec<UnsignedInt *> cell_pkg_index_set_;
     StdVec<std::pair<Arrayi, int> *> pkg_cell_info_set_;
-    StdVec<Real> global_h_ratio_vec_;
+    StdVec<Real> global_h_ratio_vec_; /**< the ratio of the reference spacing to the data spacing */
+    StdVec<NeighborMethod<SingleValued> *> neighbor_method_set_;
     StdVec<MeshWithGridDataPackagesType *> mesh_data_set_;
     StdVec<ProbeSignedDistance *> probe_signed_distance_set_;
     StdVec<ProbeNormalDirection *> probe_normal_direction_set_;
@@ -147,13 +148,13 @@ class MultilevelLevelSet : public BaseMeshField
     UniquePtrsKeeper<ProbeKernelGradientIntegral> probe_kernel_gradient_integral_vector_keeper_;
     UniquePtrsKeeper<ProbeKernelSecondGradientIntegral> probe_kernel_second_gradient_integral_vector_keeper_;
 
-    UniquePtr<BaseExecDynamics> correct_topology_keeper_;
-    UniquePtr<BaseExecDynamics> clean_interface_keeper_;
-    UniquePtr<SingularVariable<KernelTabulatedCK>> kernel_;
+    UniquePtr<BaseDynamics<void>> correct_topology_keeper_;
+    UniquePtr<BaseDynamics<void>> clean_interface_keeper_;
+    UniquePtrsKeeper<NeighborMethod<SingleValued>> neighbor_method_keeper_;
     std::function<void()> sync_mesh_variable_data_;
 
     template <class ExecutionPolicy>
-    void configLevelSetPostProcesses(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel);
+    void configLevelSetPostProcesses(const ExecutionPolicy &ex_policy);
 };
 } // namespace SPH
 #endif // LEVEL_SET_H

src/shared/geometries/level_set.hpp

@@ -7,16 +7,16 @@ namespace SPH
 {
 //=================================================================================================//
 template <class ExecutionPolicy>
-void MultilevelLevelSet::configLevelSetPostProcesses(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel)
+void MultilevelLevelSet::configLevelSetPostProcesses(const ExecutionPolicy &ex_policy)
 {
     clean_interface_keeper_ = makeUnique<CleanInterface<ExecutionPolicy>>(
-        *mesh_data_set_.back(), kernel, global_h_ratio_vec_.back());
+        *mesh_data_set_.back(), *neighbor_method_set_.back());
     correct_topology_keeper_ = makeUnique<CorrectTopology<ExecutionPolicy>>(
-        *mesh_data_set_.back(), kernel, global_h_ratio_vec_.back());
+        *mesh_data_set_.back(), *neighbor_method_set_.back());
 }
 //=================================================================================================//
 template <class ExecutionPolicy>
-void MultilevelLevelSet::initializeMeshVariables(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel)
+void MultilevelLevelSet::initializeMeshVariables(const ExecutionPolicy &ex_policy)
 {
     for (size_t level = 0; level < total_levels_; level++)
     {
@@ -35,23 +35,23 @@ void MultilevelLevelSet::initializeMeshVariables(const ExecutionPolicy &ex_polic
 template <class ExecutionPolicy>
 void MultilevelLevelSet::finishInitialization(const ExecutionPolicy &ex_policy, UsageType usage_type)
 {
-    initializeMeshVariables(ex_policy, kernel_->Data());
+    initializeMeshVariables(ex_policy);
     if (usage_type == UsageType::Volumetric)
     {
-        initializeKernelIntegralVariables(ex_policy, kernel_->Data());
-        configLevelSetPostProcesses(ex_policy, kernel_->Data());
+        initializeKernelIntegralVariables(ex_policy);
+        configLevelSetPostProcesses(ex_policy);
     }
     sync_mesh_variable_data_ = [&]()
     { this->syncMeshVariableData(ex_policy); };
 }
 //=================================================================================================//
 template <class ExecutionPolicy>
-void MultilevelLevelSet::initializeKernelIntegralVariables(const ExecutionPolicy &ex_policy, KernelTabulatedCK *kernel)
+void MultilevelLevelSet::initializeKernelIntegralVariables(const ExecutionPolicy &ex_policy)
 {
     for (size_t level = 0; level < total_levels_; level++)
     {
         MeshInnerDynamics<ExecutionPolicy, UpdateKernelIntegrals>
-            update_kernel_integrals{*mesh_data_set_[level], kernel, global_h_ratio_vec_[level]};
+            update_kernel_integrals{*mesh_data_set_[level], *neighbor_method_set_[level]};
         update_kernel_integrals.exec();
 
         probe_kernel_integral_set_.push_back(