Merge pull request #566 from Xiangyu-Hu/refactory/revise_limiter_functor

using common functors for which works for not only particles

src/shared/common/base_data_package.h

@@ -29,6 +29,7 @@
 #define BASE_DATA_PACKAGE_H
 
 #include "array_allocation.h"
+#include "common_functors.h"
 #include "data_type.h"
 #include "large_data_containers.h"
 #include "ownership.h"
@@ -80,15 +81,15 @@ class DataAssembleOperation
 
   public:
     template <typename... Args>
-    DataAssembleOperation(Args &&...args)
+    DataAssembleOperation(Args &&... args)
         : scalar_operation(std::forward<Args>(args)...),
           vector2d_operation(std::forward<Args>(args)...),
           vector3d_operation(std::forward<Args>(args)...),
           matrix2d_operation(std::forward<Args>(args)...),
           matrix3d_operation(std::forward<Args>(args)...),
           integer_operation(std::forward<Args>(args)...){};
     template <typename... OperationArgs>
-    void operator()(OperationArgs &&...operation_args)
+    void operator()(OperationArgs &&... operation_args)
     {
         scalar_operation(std::forward<OperationArgs>(operation_args)...);
         vector2d_operation(std::forward<OperationArgs>(operation_args)...);
@@ -108,18 +109,18 @@ class OperationOnDataAssemble
     OperationType operation_;
 
     template <std::size_t... Is, typename... OperationArgs>
-    void operationSequence(std::index_sequence<Is...>, OperationArgs &&...operation_args)
+    void operationSequence(std::index_sequence<Is...>, OperationArgs &&... operation_args)
     {
         (operation_(std::get<Is>(data_assemble_), std::forward<OperationArgs>(operation_args)...), ...);
     }
 
   public:
     template <typename... Args>
-    OperationOnDataAssemble(DataAssembleType &data_assemble, Args &&...args)
+    OperationOnDataAssemble(DataAssembleType &data_assemble, Args &&... args)
         : data_assemble_(data_assemble), operation_(std::forward<Args>(args)...){};
 
     template <typename... OperationArgs>
-    void operator()(OperationArgs &&...operation_args)
+    void operator()(OperationArgs &&... operation_args)
     {
         operationSequence(std::make_index_sequence<tuple_size_>{}, std::forward<OperationArgs>(operation_args)...);
     }

src/shared/common/common_functors.h

@@ -0,0 +1,80 @@
+/* ------------------------------------------------------------------------- *
+ *                                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-2023 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    common_functors.h
+ * @brief 	TBD.
+ * @author	Xiangyu Hu
+ */
+
+#ifndef COMMON_FUNCTORS_H
+#define COMMON_FUNCTORS_H
+
+#include "base_data_type.h"
+#include "scalar_functions.h"
+
+namespace SPH
+{
+
+/** 
+ * @class Limiter
+ * Base class introduce the concept of limiter, 
+ * which limits the magnitude of a value with a fraction. 
+ * The derived class should implement the operator Real() 
+ * to indicate limiting fraction.
+ * Generally, the object of the derived class 
+ * should be named as "limiter" or "limiter_" (class member) 
+ * so that the code can be more readable.
+ */
+class Limiter
+{
+};
+
+class NoLimiter : public Limiter
+{
+  public:
+    template <typename... Args>
+    NoLimiter(Args &&... args) : Limiter(){};
+
+    template <typename... Args>
+    Real operator()(Args &&... args)
+    {
+        return 1.0;
+    };
+};
+
+class TruncatedLinear : public Limiter
+{
+    Real ref_, slope_;
+
+  public:
+    TruncatedLinear(Real ref, Real slope = 100.0)
+        : Limiter(), ref_(ref), slope_(slope){};
+    Real operator()(Real measure)
+    {
+        Real measure_scale = measure * ref_;
+        return SMIN(slope_ * measure_scale, Real(1));
+    };
+};
+
+} // namespace SPH
+#endif // COMMON_FUNCTORS_H

src/shared/materials/riemann_solver.cpp

@@ -28,36 +28,4 @@ FluidStateOut NoRiemannSolver::
     return FluidStateOut(rho_star, v_star, p_star);
 }
 //=================================================================================================//
-FluidStateOut AcousticRiemannSolver::
-    InterfaceState(const FluidStateIn &state_i, const FluidStateIn &state_j, const Vecd &e_ij)
-{
-    FluidStateOut average_state = NoRiemannSolver::InterfaceState(state_i, state_j, e_ij);
-
-    Real ul = -e_ij.dot(state_i.vel_);
-    Real ur = -e_ij.dot(state_j.vel_);
-    Real u_jump = ul - ur;
-    Real limited_mach_number = SMIN(limiter_coeff_ * SMAX(u_jump * inv_c_ave_, Real(0)), Real(1));
-
-    Real p_star = average_state.p_ + 0.5 * rho0c0_geo_ave_ * u_jump * limited_mach_number;
-    Real u_dissipative = 0.5 * (state_i.p_ - state_j.p_) * inv_rho0c0_ave_ * limited_mach_number * limited_mach_number;
-    Vecd vel_star = average_state.vel_ - e_ij * u_dissipative;
-
-    return FluidStateOut(average_state.rho_, vel_star, p_star);
-}
-//=================================================================================================//
-Real AcousticRiemannSolver::DissipativePJump(const Real &u_jump)
-{
-    return rho0c0_geo_ave_ * u_jump * SMIN(limiter_coeff_ * SMAX(u_jump * inv_c_ave_, Real(0)), Real(1));
-}
-//=================================================================================================//
-Real AcousticRiemannSolver::DissipativeUJump(const Real &p_jump)
-{
-    return p_jump * inv_rho0c0_ave_;
-}
-//=================================================================================================//
-Real DissipativeRiemannSolver::DissipativePJump(const Real &u_jump)
-{
-    return rho0c0_geo_ave_ * u_jump;
-}
-//=================================================================================================//
 } // namespace SPH

src/shared/materials/riemann_solver.h

@@ -79,34 +79,48 @@ class NoRiemannSolver
     Real rho0c0_i_, rho0c0_j_, inv_rho0c0_sum_;
 };
 
-class AcousticRiemannSolver : public NoRiemannSolver
+template <typename LimiterType>
+class BaseAcousticRiemannSolver : public NoRiemannSolver
 {
   public:
     template <class FluidI, class FluidJ>
-    AcousticRiemannSolver(FluidI &fluid_i, FluidJ &fluid_j, const Real limiter_coeff = 3.0)
+    BaseAcousticRiemannSolver(FluidI &fluid_i, FluidJ &fluid_j, Real limiter_coeff = 3.0)
         : NoRiemannSolver(fluid_i, fluid_j),
           inv_rho0c0_ave_(2.0 * inv_rho0c0_sum_),
           rho0c0_geo_ave_(2.0 * rho0c0_i_ * rho0c0_j_ * inv_rho0c0_sum_),
-          inv_c_ave_(0.5 * (rho0_i_ + rho0_j_) * inv_rho0c0_ave_),
-          limiter_coeff_(limiter_coeff){};
-    Real DissipativePJump(const Real &u_jump);
-    Real DissipativeUJump(const Real &p_jump);
-    FluidStateOut InterfaceState(const FluidStateIn &state_i, const FluidStateIn &state_j, const Vecd &e_ij);
+          limiter_(0.5 * (rho0_i_ + rho0_j_) * inv_rho0c0_ave_, limiter_coeff){};
+    Real DissipativePJump(const Real &u_jump)
+    {
+        return rho0c0_geo_ave_ * u_jump * limiter_(SMAX(u_jump, Real(0)));
+    };
+    Real DissipativeUJump(const Real &p_jump)
+    {
+        return p_jump * inv_rho0c0_ave_;
+    };
+
+    FluidStateOut InterfaceState(const FluidStateIn &state_i, const FluidStateIn &state_j, const Vecd &e_ij)
+    {
+        FluidStateOut average_state = NoRiemannSolver::InterfaceState(state_i, state_j, e_ij);
+
+        Real ul = -e_ij.dot(state_i.vel_);
+        Real ur = -e_ij.dot(state_j.vel_);
+        Real u_jump = ul - ur;
+        Real limited_mach_number = limiter_(SMAX(u_jump, Real(0)));
+
+        Real p_star = average_state.p_ + 0.5 * rho0c0_geo_ave_ * u_jump * limited_mach_number;
+        Real u_dissipative = 0.5 * (state_i.p_ - state_j.p_) * inv_rho0c0_ave_ * limited_mach_number * limited_mach_number;
+        Vecd vel_star = average_state.vel_ - e_ij * u_dissipative;
+
+        return FluidStateOut(average_state.rho_, vel_star, p_star);
+    };
 
   protected:
     Real inv_rho0c0_ave_, rho0c0_geo_ave_;
-    Real inv_c_ave_;
+    LimiterType limiter_;
     Real limiter_coeff_;
 };
+using AcousticRiemannSolver = BaseAcousticRiemannSolver<TruncatedLinear>;
+using DissipativeRiemannSolver = BaseAcousticRiemannSolver<NoLimiter>;
 
-class DissipativeRiemannSolver : public AcousticRiemannSolver
-{
-  public:
-    template <class FluidI, class FluidJ>
-    DissipativeRiemannSolver(FluidI &fluid_i, FluidJ &fluid_j)
-        : AcousticRiemannSolver(fluid_i, fluid_j){};
-    Real DissipativePJump(const Real &u_jump);
-};
 } // namespace SPH
-
 #endif // RIEMANN_SOLVER_H

src/shared/particle_dynamics/all_particle_dynamics.h

@@ -30,5 +30,5 @@
 #define ALL_PARTICLE_DYNAMICS_H
 
 #include "particle_dynamics_algorithms.h"
-
+#include "particle_functors.h"
 #endif // ALL_PARTICLE_DYNAMICS_H

src/shared/particle_dynamics/base_local_dynamics.h

@@ -32,7 +32,6 @@
 
 #include "base_data_package.h"
 #include "base_particle_dynamics.h"
-#include "particle_functors.h"
 #include "sph_data_containers.h"
 
 namespace SPH
@@ -108,7 +107,7 @@ class Average : public ReduceSumType
 {
   public:
     template <class DynamicsIdentifier, typename... Args>
-    Average(DynamicsIdentifier &identifier, Args &&...args)
+    Average(DynamicsIdentifier &identifier, Args &&... args)
         : ReduceSumType(identifier, std::forward<Args>(args)...){};
     virtual ~Average(){};
     using ReturnType = typename ReduceSumType::ReturnType;
@@ -132,7 +131,7 @@ struct ConstructorArgs
     BodyRelationType &body_relation_;
     std::tuple<OtherArgs...> others_;
     SPHBody &getSPHBody() { return body_relation_.getSPHBody(); };
-    ConstructorArgs(BodyRelationType &body_relation, OtherArgs &&...other_args)
+    ConstructorArgs(BodyRelationType &body_relation, OtherArgs &&... other_args)
         : body_relation_(body_relation), others_(std::forward<OtherArgs>(other_args)...){};
 };
 
@@ -165,7 +164,7 @@ class ComplexInteraction<LocalDynamicsName<FirstInteraction, OtherInteractions..
   public:
     template <class FirstParameterSet, typename... OtherParameterSets>
     explicit ComplexInteraction(FirstParameterSet &&first_parameter_set,
-                                OtherParameterSets &&...other_parameter_sets)
+                                OtherParameterSets &&... other_parameter_sets)
         : LocalDynamicsName<FirstInteraction, CommonParameters...>(first_parameter_set),
           other_interactions_(std::forward<OtherParameterSets>(other_parameter_sets)...){};
 

src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h

@@ -57,7 +57,7 @@ class TransportVelocityCorrection<Base, DataDelegationType, KernelCorrectionType
   protected:
     StdLargeVec<Vecd> &zero_gradient_residue_;
     KernelCorrectionType kernel_correction_;
-    ParticleScope checkWithinScope;
+    ParticleScope within_scope_;
 };
 
 template <class ResolutionType, class LimiterType, typename... CommonControlTypes>
@@ -75,7 +75,7 @@ class TransportVelocityCorrection<Inner<ResolutionType, LimiterType>, CommonCont
 
   protected:
     const Real h_ref_, correction_scaling_;
-    StdLargeVec<Real>& Vol_;
+    StdLargeVec<Real> &Vol_;
     StdLargeVec<Vecd> &pos_;
     ResolutionType h_ratio_;
     LimiterType limiter_;
@@ -93,8 +93,8 @@ class TransportVelocityCorrection<Contact<Boundary>, CommonControlTypes...>
     virtual ~TransportVelocityCorrection(){};
     void interaction(size_t index_i, Real dt = 0.0);
 
-protected:
-    StdVec<StdLargeVec<Real>*> wall_Vol_;
+  protected:
+    StdVec<StdLargeVec<Real> *> wall_Vol_;
 };
 
 template <class KernelCorrectionType, typename... CommonControlTypes>
@@ -108,7 +108,7 @@ class TransportVelocityCorrection<Contact<>, KernelCorrectionType, CommonControl
 
   protected:
     StdVec<KernelCorrectionType> contact_kernel_corrections_;
-    StdVec<StdLargeVec<Real>*> contact_Vol_;
+    StdVec<StdLargeVec<Real> *> contact_Vol_;
 };
 
 template <class ResolutionType, class LimiterType, typename... CommonControlTypes>
@@ -121,7 +121,7 @@ using TransportVelocityCorrectionComplex =
 
 template <class ParticleScope>
 using TransportVelocityLimitedCorrectionComplex =
-    BaseTransportVelocityCorrectionComplex<SingleResolution, ZeroGradientLimiter, NoKernelCorrection, ParticleScope>;
+    BaseTransportVelocityCorrectionComplex<SingleResolution, TruncatedLinear, NoKernelCorrection, ParticleScope>;
 
 template <class ParticleScope>
 using TransportVelocityCorrectionComplexAdaptive =