Rkcg for open boundary flow pull request

src/shared/particle_dynamics/fluid_dynamics/density_summation.hpp

@@ -43,6 +43,8 @@ void DensitySummation<Inner<NearSurfaceType, SummationType...>>::update(size_t i
         isNearFreeSurface(index_i)
             ? near_surface_rho_(this->rho_sum_[index_i], this->rho0_, this->rho_[index_i])
             : this->rho_sum_[index_i];
+    //** Temporary treatment *
+    this->Vol_[index_i] = this->mass_[index_i] / this->rho_[index_i];
 }
 //=================================================================================================//
 template <typename NearSurfaceType, typename... SummationType>

src/shared/particle_dynamics/fluid_dynamics/fluid_integration.h

@@ -143,6 +143,7 @@ using Integration1stHalfWithWall = ComplexInteraction<Integration1stHalf<Inner<>
 using Integration1stHalfWithWallNoRiemann = Integration1stHalfWithWall<NoRiemannSolver, NoKernelCorrection>;
 using Integration1stHalfWithWallRiemann = Integration1stHalfWithWall<AcousticRiemannSolver, NoKernelCorrection>;
 using Integration1stHalfCorrectionWithWallRiemann = Integration1stHalfWithWall<AcousticRiemannSolver, LinearGradientCorrection>;
+using Integration1stHalfCorrectionForOpenBoundaryFlowWithWallRiemann = Integration1stHalfWithWall<AcousticRiemannSolver, LinearGradientCorrectionWithBulkScope>;
 
 using MultiPhaseIntegration1stHalfWithWallRiemann =
     ComplexInteraction<Integration1stHalf<Inner<>, Contact<>, Contact<Wall>>, AcousticRiemannSolver, NoKernelCorrection>;

src/shared/particle_dynamics/fluid_dynamics/fluid_integration.hpp

@@ -84,7 +84,7 @@ void Integration1stHalf<Inner<>, RiemannSolverType, KernelCorrectionType>::inter
         Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
         const Vecd &e_ij = inner_neighborhood.e_ij_[n];
 
-        force -= (p_[index_i] * correction_(index_j) + p_[index_j] * correction_(index_i)) * dW_ijV_j * e_ij;
+        force -= (p_[index_i] * correction_(index_j, index_i) + p_[index_j] * correction_(index_i)) * dW_ijV_j * e_ij;
         rho_dissipation += riemann_solver_.DissipativeUJump(p_[index_i] - p_[index_j]) * dW_ijV_j;
     }
     force_[index_i] += force * Vol_[index_i];

src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h

@@ -123,6 +123,14 @@ template <class ParticleScope>
 using TransportVelocityCorrectionCorrectedComplex =
     BaseTransportVelocityCorrectionComplex<SingleResolution, NoLimiter, LinearGradientCorrection, ParticleScope>;
 
+template <class ParticleScope>
+using TransportVelocityCorrectionCorrectedForOpenBoundaryFlowComplex =
+    BaseTransportVelocityCorrectionComplex<SingleResolution, NoLimiter, LinearGradientCorrectionWithBulkScope, ParticleScope>;
+
+template <class ParticleScope>
+using TransportVelocityLimitedCorrectionCorrectedForOpenBoundaryFlowComplex =
+    BaseTransportVelocityCorrectionComplex<SingleResolution, TruncatedLinear, LinearGradientCorrectionWithBulkScope, ParticleScope>;
+
 template <class ParticleScope>
 using TransportVelocityLimitedCorrectionComplex =
     BaseTransportVelocityCorrectionComplex<SingleResolution, TruncatedLinear, NoKernelCorrection, ParticleScope>;

src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.hpp

@@ -45,7 +45,7 @@ void TransportVelocityCorrection<Inner<ResolutionType, LimiterType>, CommonContr
         {
             size_t index_j = inner_neighborhood.j_[n];
             // acceleration for transport velocity
-            inconsistency -= (this->kernel_correction_(index_i) + this->kernel_correction_(index_j)) *
+            inconsistency -= (this->kernel_correction_(index_i) + this->kernel_correction_(index_j, index_i)) *
                              inner_neighborhood.dW_ij_[n] * this->Vol_[index_j] * inner_neighborhood.e_ij_[n];
         }
         this->zero_gradient_residue_[index_i] = inconsistency;

src/shared/particle_dynamics/particle_functors.h

@@ -186,6 +186,11 @@ class NoKernelCorrection : public KernelCorrection
 {
   public:
     NoKernelCorrection(BaseParticles *particles) : KernelCorrection(){};
+    Real operator()(size_t index_i, size_t index_j)
+    {
+        return 1.0;
+    };
+
     Real operator()(size_t index_i)
     {
         return 1.0;
@@ -199,13 +204,41 @@ class LinearGradientCorrection : public KernelCorrection
         : KernelCorrection(),
           B_(*particles->getVariableDataByName<Matd>("LinearGradientCorrectionMatrix")){};
 
+    Matd operator()(size_t index_i, size_t index_j)
+    {
+        return B_[index_i];
+    };
+
+    Matd operator()(size_t index_i)
+    {
+        return B_[index_i];
+    };
+
+  protected:
+    StdLargeVec<Matd> &B_;
+};
+
+class LinearGradientCorrectionWithBulkScope : public KernelCorrection
+{
+  public:
+    LinearGradientCorrectionWithBulkScope(BaseParticles *particles)
+        : KernelCorrection(),
+          B_(*particles->getVariableDataByName<Matd>("LinearGradientCorrectionMatrix")),
+          within_scope_(particles){};
+
+    Matd operator()(size_t index_j, size_t index_i)
+    {
+        return within_scope_(index_j) ? B_[index_j] : B_[index_i];
+    };
+
     Matd operator()(size_t index_i)
     {
         return B_[index_i];
     };
 
   protected:
     StdLargeVec<Matd> &B_;
+    BulkParticles within_scope_;
 };
 
 class SingleResolution

tests/extra_source_and_tests/extra_src/shared/pressure_boundary/density_correciton.h

@@ -48,7 +48,7 @@ class DensitySummationPressure<Base, DataDelegationType>
     virtual ~DensitySummationPressure(){};
 
   protected:
-    StdLargeVec<Real> &rho_, &mass_, &rho_sum_;
+    StdLargeVec<Real> &rho_, &mass_, &rho_sum_, &Vol_;
     Real rho0_, inv_sigma0_, W0_;
 };
 
@@ -77,6 +77,8 @@ class DensitySummationPressure<Inner<>> : public DensitySummationPressure<Inner<
     {
         if (buffer_particle_indicator_[index_i] == 0)
             assignDensity(index_i);
+        //** Temporary treatment *
+        this->Vol_[index_i] = this->mass_[index_i] / this->rho_[index_i];
     };
 
   protected:

tests/extra_source_and_tests/extra_src/shared/pressure_boundary/density_correciton.hpp

@@ -13,6 +13,7 @@ DensitySummationPressure<Base, DataDelegationType>::DensitySummationPressure(Bas
       rho_(*this->particles_->template getVariableDataByName<Real>("Density")),
       mass_(*this->particles_->template getVariableDataByName<Real>("Mass")),
       rho_sum_(*this->particles_->template registerSharedVariable<Real>("DensitySummation")),
+      Vol_(*this->particles_->template getVariableDataByName<Real>("VolumetricMeasure")),
       rho0_(this->sph_body_.base_material_->ReferenceDensity()),
       inv_sigma0_(1.0 / this->sph_body_.sph_adaptation_->LatticeNumberDensity()),
       W0_(this->sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd)) {}

tests/extra_source_and_tests/extra_src/shared/pressure_boundary/pressure_boundary.h

@@ -36,24 +36,26 @@ namespace SPH
 {
 namespace fluid_dynamics
 {
-template <typename TargetPressure>
-class PressureCondition : public BaseFlowBoundaryCondition
+template <typename TargetPressure, class KernelCorrectionType>
+class PressureBoundaryCondition : public BaseFlowBoundaryCondition
 {
   public:
     /** default parameter indicates prescribe pressure */
-    explicit PressureCondition(BodyAlignedBoxByCell &aligned_box_part)
+    explicit PressureBoundaryCondition(BodyAlignedBoxByCell &aligned_box_part)
         : BaseFlowBoundaryCondition(aligned_box_part),
           aligned_box_(aligned_box_part.getAlignedBoxShape()),
           alignment_axis_(aligned_box_.AlignmentAxis()),
           transform_(aligned_box_.getTransform()),
           target_pressure_(*this),
-          kernel_sum_(*particles_->getVariableDataByName<Vecd>("KernelSummation")){};
-    virtual ~PressureCondition(){};
+          kernel_sum_(*particles_->getVariableDataByName<Vecd>("KernelSummation")),
+          kernel_correction_(this->particles_){};
+    virtual ~PressureBoundaryCondition(){};
     AlignedBoxShape &getAlignedBox() { return aligned_box_; };
 
     void update(size_t index_i, Real dt = 0.0)
     {
-        vel_[index_i] += 2.0 * kernel_sum_[index_i] * target_pressure_(p_[index_i]) / rho_[index_i] * dt;
+        //vel_[index_i] += 2.0 * kernel_sum_[index_i] * target_pressure_(p_[index_i]) / rho_[index_i] * dt;
+        vel_[index_i] += 2.0 * kernel_correction_(index_i) * kernel_sum_[index_i] * target_pressure_(p_[index_i]) / rho_[index_i] * dt;
 
         Vecd frame_velocity = Vecd::Zero();
         frame_velocity[alignment_axis_] = transform_.xformBaseVecToFrame(vel_[index_i])[alignment_axis_];
@@ -66,7 +68,15 @@ class PressureCondition : public BaseFlowBoundaryCondition
     Transform &transform_;
     TargetPressure target_pressure_;
     StdLargeVec<Vecd> &kernel_sum_;
+    KernelCorrectionType kernel_correction_;
 };
+
+template <typename TargetPressure>
+using PressureCondition = PressureBoundaryCondition<TargetPressure, NoKernelCorrection>;
+
+template <typename TargetPressure>
+using PressureConditionCorrection = PressureBoundaryCondition<TargetPressure, LinearGradientCorrection>;
+
 } // namespace fluid_dynamics
 } // namespace SPH
 #endif // PRESSURE_BOUNDARY_H

tests/extra_source_and_tests/test_2d_mixed_poiseuille_flow/mixed_poiseuille_flow.cpp

@@ -187,7 +187,7 @@ int main(int ac, char *av[])
     SimpleDynamics<NormalDirectionFromBodyShape> wall_boundary_normal_direction(wall_boundary);
     InteractionDynamics<NablaWVComplex> kernel_summation(water_block_inner, water_block_contact);
     InteractionWithUpdate<SpatialTemporalFreeSurfaceIndicationComplex> boundary_indicator(water_block_inner, water_block_contact);
-
+    
     Dynamics1Level<fluid_dynamics::Integration1stHalfWithWallRiemann> pressure_relaxation(water_block_inner, water_block_contact);
     Dynamics1Level<fluid_dynamics::Integration2ndHalfWithWallRiemann> density_relaxation(water_block_inner, water_block_contact);
     InteractionWithUpdate<fluid_dynamics::ViscousForceWithWall> viscous_acceleration(water_block_inner, water_block_contact);
@@ -206,8 +206,10 @@ int main(int ac, char *av[])
     fluid_dynamics::BidirectionalBuffer<RightInflowPressure> right_emitter_inflow_injection(right_emitter, in_outlet_particle_buffer);
 
     InteractionWithUpdate<fluid_dynamics::DensitySummationPressureComplex> update_fluid_density(water_block_inner, water_block_contact);
+
     SimpleDynamics<fluid_dynamics::PressureCondition<LeftInflowPressure>> left_inflow_pressure_condition(left_emitter);
     SimpleDynamics<fluid_dynamics::PressureCondition<RightInflowPressure>> right_inflow_pressure_condition(right_emitter);
+
     SimpleDynamics<fluid_dynamics::InflowVelocityCondition<InflowVelocity>> inflow_velocity_condition(left_disposer);
     //----------------------------------------------------------------------
     //	Define the methods for I/O operations, observations
@@ -264,6 +266,7 @@ int main(int ac, char *av[])
             time_instance = TickCount::now();
             Real Dt = get_fluid_advection_time_step_size.exec();
             update_fluid_density.exec();
+
             viscous_acceleration.exec();
             transport_velocity_correction.exec();
             interval_computing_time_step += TickCount::now() - time_instance;