diff --git a/src/shared/particle_dynamics/base_local_dynamics.h b/src/shared/particle_dynamics/base_local_dynamics.h index 743293720c..152e1e7a0a 100644 --- a/src/shared/particle_dynamics/base_local_dynamics.h +++ b/src/shared/particle_dynamics/base_local_dynamics.h @@ -36,33 +36,76 @@ namespace SPH { -/** A Functor for Summation */ +//---------------------------------------------------------------------- +// Particle group scope functors +//---------------------------------------------------------------------- +class AllParticles +{ + public: + AllParticles(BaseParticles *base_particles){}; + bool operator()(size_t index_i) + { + return true; + }; +}; + +template +class IndicatedParticles +{ + StdLargeVec &indicator_; + + public: + IndicatedParticles(BaseParticles *base_particles) + : indicator_(*base_particles->getVariableByName("Indicator")){}; + bool operator()(size_t index_i) + { + return indicator_[index_i] == INDICATOR; + }; +}; + +template +class NotIndicatedParticles +{ + StdLargeVec &indicator_; + + public: + NotIndicatedParticles(BaseParticles *base_particles) + : indicator_(*base_particles->getVariableByName("Indicator")){}; + bool operator()(size_t index_i) + { + return indicator_[index_i] != INDICATOR; + }; +}; + +//---------------------------------------------------------------------- +// Particle reduce functors +//---------------------------------------------------------------------- template struct ReduceSum { ReturnType operator()(const ReturnType &x, const ReturnType &y) const { return x + y; }; }; -/** A Functor for Maximum */ + struct ReduceMax { Real operator()(Real x, Real y) const { return SMAX(x, y); }; }; -/** A Functor for Minimum */ + struct ReduceMin { Real operator()(Real x, Real y) const { return SMIN(x, y); }; }; -/** A Functor for OR operator */ + struct ReduceOR { bool operator()(bool x, bool y) const { return x || y; }; }; -/** A Functor for AND operator */ + struct ReduceAND { bool operator()(bool x, bool y) const { return x && y; }; }; -/** A Functor for lower bound */ + struct ReduceLowerBound { Vecd operator()(const Vecd &x, const Vecd &y) const @@ -73,7 +116,6 @@ struct ReduceLowerBound return lower_bound; }; }; -/** A Functor for upper bound */ struct ReduceUpperBound { Vecd operator()(const Vecd &x, const Vecd &y) const diff --git a/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h b/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h index 795e63035e..7da0f6a7ce 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h +++ b/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h @@ -40,3 +40,4 @@ #include "fluid_dynamics_multi_phase.hpp" #include "fluid_surface_complex.h" #include "fluid_surface_inner.hpp" +#include "transport_velocity_correction.h" diff --git a/src/shared/particle_dynamics/fluid_dynamics/base_fluid_dynamics.h b/src/shared/particle_dynamics/fluid_dynamics/base_fluid_dynamics.h new file mode 100644 index 0000000000..6e0ae9ae75 --- /dev/null +++ b/src/shared/particle_dynamics/fluid_dynamics/base_fluid_dynamics.h @@ -0,0 +1,47 @@ +/* ------------------------------------------------------------------------- * + * 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 base_fluid_dynamics.h + * @brief Collection of headers and types used by all fluid dynamics classes. + * @author Xiangyu Hu + */ + +#ifndef BASE_FLUID_DYNAMICS_H +#define BASE_FLUID_DYNAMICS_H + +#include "all_body_relations.h" +#include "all_particle_dynamics.h" +#include "base_particles.hpp" +#include "fluid_body.h" + +namespace SPH +{ +namespace fluid_dynamics +{ +typedef DataDelegateSimple FluidDataSimple; +typedef DataDelegateInner FluidDataInner; +typedef DataDelegateContact FluidContactData; +typedef DataDelegateContact FluidContactOnly; +} // namespace fluid_dynamics +} // namespace SPH +#endif // BASE_FLUID_DYNAMICS_H diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_boundary.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_boundary.h index b27d607eca..da6f00c7cb 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_boundary.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_boundary.h @@ -126,7 +126,7 @@ class FreeStreamVelocityCorrection : public LocalDynamics, public FluidDataSimpl Real rho0_; StdLargeVec &rho_sum_; StdLargeVec &pos_, &vel_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; TargetVelocity target_velocity; public: @@ -134,13 +134,13 @@ class FreeStreamVelocityCorrection : public LocalDynamics, public FluidDataSimpl : LocalDynamics(sph_body), FluidDataSimple(sph_body), transform_(transform), rho0_(DynamicCast(this, particles_->getBaseMaterial()).ReferenceDensity()), rho_sum_(*particles_->getVariableByName("DensitySummation")), pos_(particles_->pos_), vel_(particles_->vel_), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + indicator_(*particles_->getVariableByName("Indicator")), target_velocity(*this){}; virtual ~FreeStreamVelocityCorrection(){}; void update(size_t index_i, Real dt = 0.0) { - if (surface_indicator_[index_i] == 1) + if (indicator_[index_i] == 1) { Vecd frame_position = transform_.shiftBaseStationToFrame(pos_[index_i]); Vecd frame_velocity = transform_.xformBaseVecToFrame(vel_[index_i]); diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.h index f8e2a2a027..0c38f7d09e 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.h @@ -31,6 +31,8 @@ #ifndef FLUID_DYNAMICS_COMPLEX_H #define FLUID_DYNAMICS_COMPLEX_H +#include "base_fluid_dynamics.h" + #include "fluid_dynamics_inner.h" #include "fluid_dynamics_inner.hpp" #include "solid_body.h" @@ -40,10 +42,7 @@ namespace SPH { namespace fluid_dynamics { -typedef DataDelegateContact - FluidWallData; -typedef DataDelegateContact - FluidContactData; +typedef DataDelegateContact FluidWallData; typedef DataDelegateContact FSIContactData; /** * @class InteractionWithWall @@ -74,7 +73,7 @@ class InteractionWithWall : public BaseIntegrationType, public FluidWallData */ template class BaseDensitySummationComplex - : public BaseInteractionComplex + : public BaseInteractionComplex { public: template @@ -138,40 +137,6 @@ class BaseViscousAccelerationWithWall : public InteractionWithWall; -/** - * @class TransportVelocityCorrectionComplex - * @brief transport velocity correction considering the contribution from contact bodies - */ -class TransportVelocityCorrectionComplex - : public BaseInteractionComplex -{ - public: - template - TransportVelocityCorrectionComplex(Args &&...args) - : BaseInteractionComplex( - std::forward(args)...){}; - virtual ~TransportVelocityCorrectionComplex(){}; - - inline void interaction(size_t index_i, Real dt = 0.0); -}; - -/** - * @class TransportVelocityCorrectionComplexAdaptive - * @brief transport velocity correction considering the contribution from contact bodies - */ -class TransportVelocityCorrectionComplexAdaptive - : public BaseInteractionComplex -{ - public: - template - TransportVelocityCorrectionComplexAdaptive(Args &&...args) - : BaseInteractionComplex( - std::forward(args)...){}; - virtual ~TransportVelocityCorrectionComplexAdaptive(){}; - - inline void interaction(size_t index_i, Real dt = 0.0); -}; - /** * @class BaseIntegration1stHalfWithWall * @brief template class pressure relaxation scheme together with wall boundary diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp index 30c0cf4566..04ce493687 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp @@ -32,55 +32,6 @@ void DensitySummationComplexAdaptive:: sph_adaptation_.NumberDensityScaleFactor(h_ratio_[index_i]); } //=================================================================================================// -void TransportVelocityCorrectionComplex:: - interaction(size_t index_i, Real dt) -{ - TransportVelocityCorrectionInner::interaction(index_i, dt); - - Vecd acceleration_trans = Vecd::Zero(); - for (size_t k = 0; k < contact_configuration_.size(); ++k) - { - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; - for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) - { - Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; - - // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; - } - } - - /** correcting particle position */ - if (surface_indicator_[index_i] == 0) - pos_[index_i] += coefficient_ * smoothing_length_sqr_ * acceleration_trans; -} -//=================================================================================================// -void TransportVelocityCorrectionComplexAdaptive:: - interaction(size_t index_i, Real dt) -{ - TransportVelocityCorrectionInnerAdaptive::interaction(index_i, dt); - - Vecd acceleration_trans = Vecd::Zero(); - for (size_t k = 0; k < contact_configuration_.size(); ++k) - { - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; - for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) - { - Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; - - // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; - } - } - - /** correcting particle position */ - if (surface_indicator_[index_i] == 0) - { - Real inv_h_ratio = 1.0 / sph_adaptation_.SmoothingLengthRatio(index_i); - pos_[index_i] += coefficient_ * smoothing_length_sqr_ * inv_h_ratio * inv_h_ratio * acceleration_trans; - } -} -//=================================================================================================// void Oldroyd_BIntegration1stHalfWithWall:: interaction(size_t index_i, Real dt) { @@ -160,7 +111,7 @@ template template BaseDensitySummationComplex:: BaseDensitySummationComplex(Args &&...args) - : BaseInteractionComplex(std::forward(args)...) + : BaseInteractionComplex(std::forward(args)...) { for (size_t k = 0; k != this->contact_particles_.size(); ++k) { diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp index 3d3acbff5b..3dac7ba91d 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp @@ -42,21 +42,6 @@ BaseViscousAccelerationInner::BaseViscousAccelerationInner(BaseInnerRelation &in mu_(DynamicCast(this, particles_->getBaseMaterial()).ReferenceViscosity()), smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()) {} //=================================================================================================// -TransportVelocityCorrectionInner:: - TransportVelocityCorrectionInner(BaseInnerRelation &inner_relation, Real coefficient) - : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), - pos_(particles_->pos_), surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), - smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)), - coefficient_(coefficient) {} -//=================================================================================================// -TransportVelocityCorrectionInnerAdaptive:: - TransportVelocityCorrectionInnerAdaptive(BaseInnerRelation &inner_relation, Real coefficient) - : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), - sph_adaptation_(*sph_body_.sph_adaptation_), - pos_(particles_->pos_), surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), - smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)), - coefficient_(coefficient) {} -//=================================================================================================// AcousticTimeStepSize::AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL) : LocalDynamicsReduce(sph_body, Real(0)), FluidDataSimple(sph_body), fluid_(DynamicCast(this, particles_->getBaseMaterial())), diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h index aa7f480b26..b4ea5c1fc9 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h @@ -32,11 +32,8 @@ #ifndef FLUID_DYNAMICS_INNER_H #define FLUID_DYNAMICS_INNER_H -#include "all_body_relations.h" -#include "all_particle_dynamics.h" -#include "base_kernel.h" -#include "base_particles.hpp" -#include "fluid_body.h" +#include "base_fluid_dynamics.h" + #include "riemann_solver.h" #include "weakly_compressible_fluid.h" @@ -44,9 +41,6 @@ namespace SPH { namespace fluid_dynamics { -typedef DataDelegateSimple FluidDataSimple; -typedef DataDelegateInner FluidDataInner; - /** * @class FluidInitialCondition * @brief Set initial condition for a fluid body. @@ -158,52 +152,6 @@ class AngularConservativeViscousAccelerationInner : public BaseViscousAccelerati inline void interaction(size_t index_i, Real dt = 0.0); }; -/** - * @class TransportVelocityCorrectionInner - * @brief The particle positions are corrected for more uniformed distribution - * when there is negative pressure in the flow. - * @details Note that the default coefficient is for using the dual time criteria method: - * Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics. - * C Zhang, M Rezavand, X Hu - Journal of Computational Physics, - * Volume 404, 1 March 2020, 109135. - * If single (acoustic) time step is used, the coefficient should be decrease - * to about 1/4 of the default value. - */ -class TransportVelocityCorrectionInner : public LocalDynamics, public FluidDataInner -{ - public: - explicit TransportVelocityCorrectionInner(BaseInnerRelation &inner_relation, Real coefficient = 0.2); - virtual ~TransportVelocityCorrectionInner(){}; - - inline void interaction(size_t index_i, Real dt = 0.0); - - protected: - StdLargeVec &pos_; - StdLargeVec &surface_indicator_; - Real smoothing_length_sqr_; - const Real coefficient_; -}; - -/** - * @class TransportVelocityCorrectionInner - * @brief transport velocity correction - */ -class TransportVelocityCorrectionInnerAdaptive : public LocalDynamics, public FluidDataInner -{ - public: - explicit TransportVelocityCorrectionInnerAdaptive(BaseInnerRelation &inner_relation, Real coefficient = 0.2); - virtual ~TransportVelocityCorrectionInnerAdaptive(){}; - - inline void interaction(size_t index_i, Real dt = 0.0); - - protected: - SPHAdaptation &sph_adaptation_; - StdLargeVec &pos_; - StdLargeVec &surface_indicator_; - Real smoothing_length_sqr_; - const Real coefficient_; -}; - /** * @class AcousticTimeStepSize * @brief Computing the acoustic time step size diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.hpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.hpp index 2a3e69f1b1..791186224c 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.hpp @@ -77,23 +77,6 @@ void AngularConservativeViscousAccelerationInner:: acc_prior_[index_i] += acceleration / rho_[index_i]; } //=================================================================================================// -void TransportVelocityCorrectionInner:: - interaction(size_t index_i, Real dt) -{ - Vecd acceleration_trans = Vecd::Zero(); - const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; - for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) - { - Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; - - // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; - } - - if (surface_indicator_[index_i] == 0) - pos_[index_i] += coefficient_ * smoothing_length_sqr_ * acceleration_trans; -} -//=================================================================================================// void VorticityInner::interaction(size_t index_i, Real dt) { AngularVecd vorticity = ZeroData::value; @@ -128,26 +111,6 @@ void Oldroyd_BIntegration1stHalf:: acc_[index_i] += acceleration / rho_[index_i]; } //=================================================================================================// -void TransportVelocityCorrectionInnerAdaptive:: - interaction(size_t index_i, Real dt) -{ - Vecd acceleration_trans = Vecd::Zero(); - const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; - for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) - { - Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; - - // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; - } - - if (surface_indicator_[index_i] == 0) - { - Real inv_h_ratio = 1.0 / sph_adaptation_.SmoothingLengthRatio(index_i); - pos_[index_i] += coefficient_ * smoothing_length_sqr_ * inv_h_ratio * inv_h_ratio * acceleration_trans; - } -} -//=================================================================================================// void Oldroyd_BIntegration2ndHalf:: interaction(size_t index_i, Real dt) { diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp index f61e73d53f..5f4c130cb1 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp @@ -35,7 +35,7 @@ MultiPhaseColorFunctionGradient:: : LocalDynamics(contact_relation.getSPHBody()), MultiPhaseData(contact_relation), rho0_(sph_body_.base_material_->ReferenceDensity()), Vol_(particles_->Vol_), pos_div_(*particles_->getVariableByName("PositionDivergence")), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")) + indicator_(*particles_->getVariableByName("Indicator")) { particles_->registerVariable(color_grad_, "ColorGradient"); particles_->registerVariable(surface_norm_, "SurfaceNormal"); diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h index 4701c401d7..1f6ffa6a81 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h @@ -149,7 +149,7 @@ class MultiPhaseColorFunctionGradient : public LocalDynamics, public MultiPhaseD Real rho0_; StdVec contact_rho0_; StdLargeVec &Vol_, &pos_div_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec color_grad_, surface_norm_; StdVec *> contact_Vol_; }; diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp index a9798af3e4..627e2a040f 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp @@ -40,7 +40,7 @@ void MultiPhaseColorFunctionGradient:: { Real Vol_i = Vol_[index_i]; Vecd gradient = Vecd::Zero(); - if (surface_indicator_[index_i]) + if (indicator_[index_i]) { for (size_t k = 0; k < contact_configuration_.size(); ++k) { diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp index d401a4f25f..20d2fced49 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp @@ -9,7 +9,7 @@ namespace fluid_dynamics FreeSurfaceIndicationComplex:: FreeSurfaceIndicationComplex(BaseInnerRelation &inner_relation, BaseContactRelation &contact_relation, Real threshold) - : FreeSurfaceIndicationInner(inner_relation, threshold), FluidContactData(contact_relation) + : FreeSurfaceIndicationInner(inner_relation, threshold), FluidContactOnly(contact_relation) { for (size_t k = 0; k != contact_particles_.size(); ++k) { @@ -26,7 +26,7 @@ FreeSurfaceIndicationComplex:: //=================================================================================================// ColorFunctionGradientComplex::ColorFunctionGradientComplex(BaseInnerRelation &inner_relation, BaseContactRelation &contact_relation) - : ColorFunctionGradientInner(inner_relation), FluidContactData(contact_relation) + : ColorFunctionGradientInner(inner_relation), FluidContactOnly(contact_relation) { for (size_t k = 0; k != contact_particles_.size(); ++k) { @@ -41,7 +41,7 @@ ColorFunctionGradientComplex::ColorFunctionGradientComplex(ComplexRelation &comp SurfaceNormWithWall::SurfaceNormWithWall(BaseContactRelation &contact_relation, Real contact_angle) : LocalDynamics(contact_relation.getSPHBody()), FSIContactData(contact_relation), contact_angle_(contact_angle), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + indicator_(*particles_->getVariableByName("Indicator")), surface_norm_(*particles_->getVariableByName("SurfaceNormal")), pos_div_(*particles_->getVariableByName("PositionDivergence")) { diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h index dedca548c3..697a1ab361 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h @@ -42,7 +42,7 @@ namespace fluid_dynamics * @class FreeSurfaceIndicationComplex * @brief indicate the particles near the free fluid surface. */ -class FreeSurfaceIndicationComplex : public FreeSurfaceIndicationInner, public FluidContactData +class FreeSurfaceIndicationComplex : public FreeSurfaceIndicationInner, public FluidContactOnly { public: FreeSurfaceIndicationComplex(BaseInnerRelation &inner_relation, @@ -84,7 +84,7 @@ using DensitySummationFreeStreamComplexAdaptive = DensitySummationFreeStream &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec &surface_norm_; StdLargeVec &pos_div_; StdVec *> wall_n_; diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.cpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.cpp index b37589ae0e..d281fd66f2 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.cpp @@ -10,7 +10,7 @@ FreeSurfaceIndicationInner:: FreeSurfaceIndicationInner(BaseInnerRelation &inner_relation, Real threshold) : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), threshold_by_dimensions_(threshold * (Real)Dimensions), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + indicator_(*particles_->getVariableByName("Indicator")), smoothing_length_(inner_relation.getSPHBody().sph_adaptation_->ReferenceSmoothingLength()) { particles_->registerVariable(pos_div_, "PositionDivergence"); @@ -18,9 +18,9 @@ FreeSurfaceIndicationInner:: //=================================================================================================// void FreeSurfaceIndicationInner::update(size_t index_i, Real dt) { - surface_indicator_[index_i] = 1; + indicator_[index_i] = 1; if (pos_div_[index_i] > threshold_by_dimensions_ && !isVeryNearFreeSurface(index_i)) - surface_indicator_[index_i] = 0; + indicator_[index_i] = 0; } //=================================================================================================// bool FreeSurfaceIndicationInner::isVeryNearFreeSurface(size_t index_i) @@ -42,7 +42,7 @@ bool FreeSurfaceIndicationInner::isVeryNearFreeSurface(size_t index_i) //=================================================================================================// ColorFunctionGradientInner::ColorFunctionGradientInner(BaseInnerRelation &inner_relation) : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + indicator_(*particles_->getVariableByName("Indicator")), pos_div_(*particles_->getVariableByName("PositionDivergence")), threshold_by_dimensions_((0.75 * (Real)Dimensions)) { @@ -52,7 +52,7 @@ ColorFunctionGradientInner::ColorFunctionGradientInner(BaseInnerRelation &inner_ //=================================================================================================// ColorFunctionGradientInterpolationInner::ColorFunctionGradientInterpolationInner(BaseInnerRelation &inner_relation) : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), Vol_(particles_->Vol_), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + indicator_(*particles_->getVariableByName("Indicator")), color_grad_(*particles_->getVariableByName("ColorGradient")), surface_norm_(*particles_->getVariableByName("SurfaceNormal")), pos_div_(*particles_->getVariableByName("PositionDivergence")), @@ -66,7 +66,7 @@ ColorFunctionGradientInterpolationInner::ColorFunctionGradientInterpolationInner SurfaceTensionAccelerationInner::SurfaceTensionAccelerationInner(BaseInnerRelation &inner_relation, Real gamma) : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), gamma_(gamma), Vol_(particles_->Vol_), mass_(particles_->mass_), - acc_prior_(particles_->acc_prior_), surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")), + acc_prior_(particles_->acc_prior_), indicator_(*particles_->getVariableByName("Indicator")), color_grad_(*particles_->getVariableByName("ColorGradient")), surface_norm_(*particles_->getVariableByName("SurfaceNormal")) {} //=================================================================================================// diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.h b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.h index 66038c2c31..cbdfa77e78 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.h +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.h @@ -55,7 +55,7 @@ class FreeSurfaceIndicationInner : public LocalDynamics, public FluidDataInner protected: Real threshold_by_dimensions_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec pos_div_; Real smoothing_length_; bool isVeryNearFreeSurface(size_t index_i); @@ -122,7 +122,7 @@ class DensitySummationFreeStream : public DensitySummationFreeSurfaceType void update(size_t index_i, Real dt = 0.0); protected: - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; bool isNearFreeSurface(size_t index_i); }; @@ -160,7 +160,7 @@ class ColorFunctionGradientInner : public LocalDynamics, public FluidDataInner inline void interaction(size_t index_i, Real dt = 0.0); protected: - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec &pos_div_; Real threshold_by_dimensions_; StdLargeVec color_grad_; @@ -181,7 +181,7 @@ class ColorFunctionGradientInterpolationInner : public LocalDynamics, public Flu protected: StdLargeVec &Vol_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec &color_grad_; StdLargeVec &surface_norm_; StdLargeVec &pos_div_; @@ -205,7 +205,7 @@ class SurfaceTensionAccelerationInner : public LocalDynamics, public FluidDataIn Real gamma_; StdLargeVec &Vol_, &mass_; StdLargeVec &acc_prior_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec &color_grad_; StdLargeVec &surface_norm_; }; diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.hpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.hpp index 749f808094..206808a149 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_inner.hpp @@ -49,13 +49,13 @@ void ColorFunctionGradientInterpolationInner:: Vecd grad = Vecd::Zero(); Real weight(0); Real total_weight(0); - if (surface_indicator_[index_i] == 1 && pos_div_[index_i] > threshold_by_dimensions_) + if (indicator_[index_i] == 1 && pos_div_[index_i] > threshold_by_dimensions_) { Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - if (surface_indicator_[index_j] == 1 && pos_div_[index_j] < threshold_by_dimensions_) + if (indicator_[index_j] == 1 && pos_div_[index_j] < threshold_by_dimensions_) { weight = inner_neighborhood.W_ij_[n] * Vol_[index_j]; grad += weight * color_grad_[index_j]; @@ -75,13 +75,13 @@ void SurfaceTensionAccelerationInner:: Real curvature(0.0); Real renormalized_curvature(0); Real pos_div(0); - if (surface_indicator_[index_i] == 1) + if (indicator_[index_i] == 1) { Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - if (surface_indicator_[index_j] == 1) + if (indicator_[index_j] == 1) { Vecd n_j = surface_norm_[index_j]; Vecd n_ij = n_i - n_j; @@ -147,7 +147,7 @@ void SpatialTemporalFreeSurfaceIdentification:: { FreeSurfaceIdentification::update(index_i, dt); - previous_surface_indicator_[index_i] = this->surface_indicator_[index_i]; + previous_surface_indicator_[index_i] = this->indicator_[index_i]; } //=================================================================================================// template @@ -161,7 +161,7 @@ template DensitySummationFreeStream:: DensitySummationFreeStream(ConstructorArgs &&...args) : DensitySummationFreeSurfaceType(std::forward(args)...), - surface_indicator_(*this->particles_->template getVariableByName("SurfaceIndicator")){}; + indicator_(*this->particles_->template getVariableByName("Indicator")){}; //=================================================================================================// template void DensitySummationFreeStream::update(size_t index_i, Real dt) @@ -183,7 +183,7 @@ bool DensitySummationFreeStream::isNearFreeSurf const Neighborhood &inner_neighborhood = this->inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { - if (surface_indicator_[inner_neighborhood.j_[n]] == 1) + if (indicator_[inner_neighborhood.j_[n]] == 1) { is_near_surface = true; break; diff --git a/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.cpp b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.cpp new file mode 100644 index 0000000000..c6fc84a2ac --- /dev/null +++ b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.cpp @@ -0,0 +1,63 @@ +#include "transport_velocity_correction.h" + +namespace SPH +{ +namespace fluid_dynamics +{ +//=================================================================================================// +TransportVelocityCorrectionInnerAdaptive:: + TransportVelocityCorrectionInnerAdaptive(BaseInnerRelation &inner_relation, Real coefficient) + : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), + sph_adaptation_(*sph_body_.sph_adaptation_), + pos_(particles_->pos_), indicator_(*particles_->getVariableByName("Indicator")), + smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)), + coefficient_(coefficient) {} + +//=================================================================================================// +void TransportVelocityCorrectionInnerAdaptive:: + interaction(size_t index_i, Real dt) +{ + Vecd acceleration_trans = Vecd::Zero(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + + // acceleration for transport velocity + acceleration_trans -= 2.0 * nablaW_ijV_j; + } + + if (indicator_[index_i] == 0) + { + Real inv_h_ratio = 1.0 / sph_adaptation_.SmoothingLengthRatio(index_i); + pos_[index_i] += coefficient_ * smoothing_length_sqr_ * inv_h_ratio * inv_h_ratio * acceleration_trans; + } +} +//=================================================================================================// +void TransportVelocityCorrectionComplexAdaptive::interaction(size_t index_i, Real dt) +{ + TransportVelocityCorrectionInnerAdaptive::interaction(index_i, dt); + + Vecd acceleration_trans = Vecd::Zero(); + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; + + // acceleration for transport velocity + acceleration_trans -= 2.0 * nablaW_ijV_j; + } + } + + /** correcting particle position */ + if (indicator_[index_i] == 0) + { + Real inv_h_ratio = 1.0 / sph_adaptation_.SmoothingLengthRatio(index_i); + pos_[index_i] += coefficient_ * smoothing_length_sqr_ * inv_h_ratio * inv_h_ratio * acceleration_trans; + } +} +//=================================================================================================// +} // namespace fluid_dynamics +} // namespace SPH diff --git a/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h new file mode 100644 index 0000000000..2dd02c3974 --- /dev/null +++ b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h @@ -0,0 +1,170 @@ +/* ------------------------------------------------------------------------- * + * 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 transport_velocity_correction.h + * @brief The particle positions are corrected for more uniformed distribution + * when there is negative pressure in the flow. + * @details Note that the default coefficient is for using the dual time criteria method: + * Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics. + * C Zhang, M Rezavand, X Hu - Journal of Computational Physics, + * Volume 404, 1 March 2020, 109135. + * If single (acoustic) time step is used, the coefficient should be decrease + * to about 1/4 of the default value. + * @author Chi Zhang and Xiangyu Hu + */ + +#ifndef TRANSPORT_VELOCITY_CORRECTION_H +#define TRANSPORT_VELOCITY_CORRECTION_H + +#include "base_fluid_dynamics.h" + +namespace SPH +{ +namespace fluid_dynamics +{ +/** + * @class TransportVelocityCorrectionInner + * @brief The particle positions are corrected for more uniformed distribution + * when there is negative pressure in the flow. + * @details Note that the default coefficient is for using the dual time criteria method: + * Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics. + * C Zhang, M Rezavand, X Hu - Journal of Computational Physics, + * Volume 404, 1 March 2020, 109135. + * If single (acoustic) time step is used, the coefficient should be decrease + * to about 1/4 of the default value. + */ +template +class TransportVelocityCorrectionInner : public LocalDynamics, public FluidDataInner +{ + public: + explicit TransportVelocityCorrectionInner(BaseInnerRelation &inner_relation, Real coefficient = 0.2) + : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation), + pos_(particles_->pos_), indicator_(*particles_->getVariableByName("Indicator")), + smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)), + coefficient_(coefficient), checkWithinScope(particles_){}; + virtual ~TransportVelocityCorrectionInner(){}; + + void interaction(size_t index_i, Real dt = 0.0) + { + if (checkWithinScope(index_i)) + { + Vecd acceleration_trans = Vecd::Zero(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + + // acceleration for transport velocity + acceleration_trans -= 2.0 * nablaW_ijV_j; + } + + pos_[index_i] += coefficient_ * smoothing_length_sqr_ * acceleration_trans; + } + }; + + protected: + StdLargeVec &pos_; + StdLargeVec &indicator_; + Real smoothing_length_sqr_; + const Real coefficient_; + ParticleScopeType checkWithinScope; +}; + +/** + * @class TransportVelocityCorrectionComplex + * @brief transport velocity correction considering the contribution from contact bodies + */ +template +class TransportVelocityCorrectionComplex + : public BaseInteractionComplex, FluidContactOnly> +{ + public: + template + TransportVelocityCorrectionComplex(Args &&...args) + : BaseInteractionComplex, FluidContactOnly>( + std::forward(args)...){}; + virtual ~TransportVelocityCorrectionComplex(){}; + + void interaction(size_t index_i, Real dt = 0.0) + { + TransportVelocityCorrectionInner::interaction(index_i, dt); + + if (this->checkWithinScope(index_i)) + { + Vecd acceleration_trans = Vecd::Zero(); + for (size_t k = 0; k < this->contact_configuration_.size(); ++k) + { + Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; + + // acceleration for transport velocity + acceleration_trans -= 2.0 * nablaW_ijV_j; + } + } + + this->pos_[index_i] += this->coefficient_ * this->smoothing_length_sqr_ * acceleration_trans; + } + }; +}; + +/** + * @class TransportVelocityCorrectionInnerAdaptive + * @brief transport velocity correction + */ +class TransportVelocityCorrectionInnerAdaptive : public LocalDynamics, public FluidDataInner +{ + public: + explicit TransportVelocityCorrectionInnerAdaptive(BaseInnerRelation &inner_relation, Real coefficient = 0.2); + virtual ~TransportVelocityCorrectionInnerAdaptive(){}; + + void interaction(size_t index_i, Real dt = 0.0); + + protected: + SPHAdaptation &sph_adaptation_; + StdLargeVec &pos_; + StdLargeVec &indicator_; + Real smoothing_length_sqr_; + const Real coefficient_; +}; + +/** + * @class TransportVelocityCorrectionComplexAdaptive + * @brief transport velocity correction considering the contribution from contact bodies + */ +class TransportVelocityCorrectionComplexAdaptive + : public BaseInteractionComplex +{ + public: + template + TransportVelocityCorrectionComplexAdaptive(Args &&...args) + : BaseInteractionComplex( + std::forward(args)...){}; + virtual ~TransportVelocityCorrectionComplexAdaptive(){}; + + void interaction(size_t index_i, Real dt = 0.0); +}; +} // namespace fluid_dynamics +} // namespace SPH +#endif // TRANSPORT_VELOCITY_CORRECTION_H diff --git a/src/shared/particles/base_particles.cpp b/src/shared/particles/base_particles.cpp index e48d521053..5ce494facb 100644 --- a/src/shared/particles/base_particles.cpp +++ b/src/shared/particles/base_particles.cpp @@ -42,7 +42,7 @@ void BaseParticles::initializeOtherVariables() registerVariable(mass_, "MassiveMeasure", [&](size_t i) -> Real { return rho_[i] * Vol_[i]; }); - registerVariable(surface_indicator_, "SurfaceIndicator"); + registerVariable(indicator_, "Indicator"); /** * add basic output particle data */ diff --git a/src/shared/particles/base_particles.h b/src/shared/particles/base_particles.h index ec9f2200e9..e7f259da25 100644 --- a/src/shared/particles/base_particles.h +++ b/src/shared/particles/base_particles.h @@ -97,10 +97,10 @@ class BaseParticles StdLargeVec acc_; /**< Acceleration induced by pressure- or stress */ StdLargeVec acc_prior_; /**< Other, such as gravity and viscous, accelerations computed before acc_ */ - StdLargeVec Vol_; /**< Volumetric measure, also area and length of surface and linear particle */ - StdLargeVec rho_; /**< Density */ - StdLargeVec mass_; /**< Massive measure, also mass per-unit thickness and per-unit cross-section of surface and linear particle */ - StdLargeVec surface_indicator_; /**< free surface indicator */ + StdLargeVec Vol_; /**< Volumetric measure, also area and length of surface and linear particle */ + StdLargeVec rho_; /**< Density */ + StdLargeVec mass_; /**< Massive measure, also mass per-unit thickness and per-unit cross-section of surface and linear particle */ + StdLargeVec indicator_; /**< particle indicator: 0 for bulk, 1 for free surface indicator, other to be defined */ //---------------------------------------------------------------------- // Global information for defining particle groups //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp index b6b5d3cdf8..6a8cabf9cd 100644 --- a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp +++ b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp @@ -125,12 +125,12 @@ int main(int ac, char *av[]) Dynamics1Level pressure_relaxation(water_block_complex_relation); Dynamics1Level density_relaxation(water_block_complex_relation); InteractionDynamics viscous_acceleration(water_block_complex_relation); - InteractionDynamics transport_velocity_correction(water_block_complex_relation); + InteractionDynamics>> transport_velocity_correction(water_block_complex_relation); InteractionWithUpdate inlet_outlet_surface_particle_indicator(water_block_complex_relation); InteractionWithUpdate update_density_by_summation(water_block_complex_relation); water_block.addBodyStateForRecording("Pressure"); // output for debug - water_block.addBodyStateForRecording("SurfaceIndicator"); // output for debug + water_block.addBodyStateForRecording("Indicator"); // output for debug SimpleDynamics initialize_a_fluid_step(water_block); ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); diff --git a/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp b/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp index c006a384b4..5f1c00d020 100644 --- a/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp +++ b/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp @@ -33,7 +33,7 @@ int main(int ac, char *av[]) (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) ? water_block.generateParticles(io_environment, water_block.getName()) : water_block.generateParticles(); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); SolidBody cylinder(sph_system, makeShared("Cylinder")); cylinder.defineAdaptationRatios(1.15, 2.0); diff --git a/tests/2d_examples/test_2d_filling_tank/filling_tank.cpp b/tests/2d_examples/test_2d_filling_tank/filling_tank.cpp index 721bd5604b..5ccb026297 100644 --- a/tests/2d_examples/test_2d_filling_tank/filling_tank.cpp +++ b/tests/2d_examples/test_2d_filling_tank/filling_tank.cpp @@ -123,7 +123,7 @@ int main() InteractionWithUpdate indicate_free_surface(water_body_complex); water_body.addBodyStateForRecording("PositionDivergence"); // for debug - water_body.addBodyStateForRecording("SurfaceIndicator"); // for debug + water_body.addBodyStateForRecording("Indicator"); // for debug SharedPtr gravity_ptr = makeShared(Vecd(0.0, -gravity_g)); SimpleDynamics initialize_a_fluid_step(water_body, gravity_ptr); diff --git a/tests/2d_examples/test_2d_flow_around_cylinder/2d_flow_around_cylinder.cpp b/tests/2d_examples/test_2d_flow_around_cylinder/2d_flow_around_cylinder.cpp index 96c1ed834b..a8b452e5a1 100644 --- a/tests/2d_examples/test_2d_flow_around_cylinder/2d_flow_around_cylinder.cpp +++ b/tests/2d_examples/test_2d_flow_around_cylinder/2d_flow_around_cylinder.cpp @@ -116,7 +116,7 @@ int main(int ac, char *av[]) /** Computing viscous acceleration with wall. */ InteractionDynamics viscous_acceleration(water_block_complex); /** Impose transport velocity. */ - InteractionDynamics transport_velocity_correction(water_block_complex); + InteractionDynamics> transport_velocity_correction(water_block_complex); /** Computing vorticity in the flow. */ InteractionDynamics compute_vorticity(water_block_complex.getInnerRelation()); /** free stream boundary condition. */ diff --git a/tests/2d_examples/test_2d_free_stream_around_cylinder/2d_free_stream_around_cylinder.cpp b/tests/2d_examples/test_2d_free_stream_around_cylinder/2d_free_stream_around_cylinder.cpp index 27f98322aa..68c03ec08f 100644 --- a/tests/2d_examples/test_2d_free_stream_around_cylinder/2d_free_stream_around_cylinder.cpp +++ b/tests/2d_examples/test_2d_free_stream_around_cylinder/2d_free_stream_around_cylinder.cpp @@ -114,7 +114,7 @@ int main(int ac, char *av[]) InteractionWithUpdate update_fluid_density(water_block_complex); /** We can output a method-specific particle data for debug */ water_block.addBodyStateForRecording("Pressure"); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); /** Time step size without considering sound wave speed. */ ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); /** Time step size with considering sound wave speed. */ @@ -130,7 +130,7 @@ int main(int ac, char *av[]) /** Computing viscous acceleration. */ InteractionDynamics viscous_acceleration(water_block_complex); /** Apply transport velocity formulation. */ - InteractionDynamics transport_velocity_correction(water_block_complex); + InteractionDynamics>> transport_velocity_correction(water_block_complex); /** compute the vorticity. */ InteractionDynamics compute_vorticity(water_block_inner); //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_free_stream_around_cylinder_mr/mr_free_stream_around_cylinder.cpp b/tests/2d_examples/test_2d_free_stream_around_cylinder_mr/mr_free_stream_around_cylinder.cpp index 62c6e9ace3..727e6e3dd1 100644 --- a/tests/2d_examples/test_2d_free_stream_around_cylinder_mr/mr_free_stream_around_cylinder.cpp +++ b/tests/2d_examples/test_2d_free_stream_around_cylinder_mr/mr_free_stream_around_cylinder.cpp @@ -127,7 +127,7 @@ int main(int ac, char *av[]) InteractionWithUpdate update_fluid_density(water_block_complex); /** We can output a method-specific particle data for debug */ water_block.addBodyStateForRecording("Pressure"); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); water_block.addBodyStateForRecording("VolumetricMeasure"); /** Time step size without considering sound wave speed. */ ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); diff --git a/tests/2d_examples/test_2d_fsi2/fsi2.cpp b/tests/2d_examples/test_2d_fsi2/fsi2.cpp index af7b994604..d1619dd8e8 100644 --- a/tests/2d_examples/test_2d_fsi2/fsi2.cpp +++ b/tests/2d_examples/test_2d_fsi2/fsi2.cpp @@ -119,7 +119,7 @@ int main(int ac, char *av[]) Dynamics1Level pressure_relaxation(water_block_complex); Dynamics1Level density_relaxation(water_block_complex); /** viscous acceleration and transport velocity correction can be combined because they are independent dynamics. */ - InteractionDynamics transport_correction(water_block_complex); + InteractionDynamics> transport_correction(water_block_complex); InteractionDynamics viscous_acceleration(water_block_complex); /** Computing vorticity in the flow. */ InteractionDynamics compute_vorticity(water_block_complex.getInnerRelation()); diff --git a/tests/2d_examples/test_2d_heat_transfer/heat_transfer.cpp b/tests/2d_examples/test_2d_heat_transfer/heat_transfer.cpp index cd46f8e423..f3fce68c46 100644 --- a/tests/2d_examples/test_2d_heat_transfer/heat_transfer.cpp +++ b/tests/2d_examples/test_2d_heat_transfer/heat_transfer.cpp @@ -274,7 +274,7 @@ int main() /** Computing viscous acceleration. */ InteractionDynamics viscous_acceleration(fluid_body_complex); /** Apply transport velocity formulation. */ - InteractionDynamics transport_velocity_correction(fluid_body_complex); + InteractionDynamics> transport_velocity_correction(fluid_body_complex); /** Computing vorticity in the flow. */ InteractionDynamics compute_vorticity(fluid_body_inner); /** Inflow boundary condition. */ diff --git a/tests/2d_examples/test_2d_poiseuille_flow/src/poiseuille_flow.cpp b/tests/2d_examples/test_2d_poiseuille_flow/src/poiseuille_flow.cpp index b40ed87371..41c851e9a7 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/src/poiseuille_flow.cpp +++ b/tests/2d_examples/test_2d_poiseuille_flow/src/poiseuille_flow.cpp @@ -127,7 +127,7 @@ int main() /** Computing viscous acceleration. */ InteractionDynamics viscous_acceleration(water_block_complex); /** Impose transport velocity. */ - InteractionDynamics transport_velocity_correction(water_block_complex); + InteractionDynamics> transport_velocity_correction(water_block_complex); /** * @brief Output. */ diff --git a/tests/2d_examples/test_2d_square_droplet/src/droplet.cpp b/tests/2d_examples/test_2d_square_droplet/src/droplet.cpp index a9bd392757..6e3d7bf78b 100644 --- a/tests/2d_examples/test_2d_square_droplet/src/droplet.cpp +++ b/tests/2d_examples/test_2d_square_droplet/src/droplet.cpp @@ -15,7 +15,7 @@ Real DH = 2.0; /**< Tank height. */ Real LL = 1.0; /**< Liquid column length. */ Real LH = 1.0; /**< Liquid column height. */ Real particle_spacing_ref = DL / 40.0; /**< Initial reference particle spacing. */ -Real BW = particle_spacing_ref * 2; /**< Extending width for BCs. */ +Real BW = particle_spacing_ref * 4; /**< Extending width for BCs. */ //---------------------------------------------------------------------- // Material parameters. //---------------------------------------------------------------------- @@ -116,7 +116,7 @@ int main() FluidBody water_block(sph_system, makeShared("WaterBody")); water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); water_block.generateParticles(); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); FluidBody air_block(sph_system, makeShared("AirBody")); air_block.defineParticlesAndMaterial(rho0_a, c_f, mu_a); @@ -146,8 +146,10 @@ int main() update_water_density_by_summation(water_wall_contact, water_air_complex); InteractionWithUpdate update_air_density_by_summation(air_wall_contact, air_water_complex); - InteractionDynamics - air_transport_correction(air_wall_contact, air_water_complex, 0.05); + InteractionDynamics> + air_transport_correction(air_wall_contact, air_water_complex, 2.0e-3); + InteractionDynamics> + water_transport_correction(water_air_complex, 2.0e-3); /** Time step size without considering sound wave speed. */ ReduceDynamics get_water_advection_time_step_size(water_block, U_ref); ReduceDynamics get_air_advection_time_step_size(air_block, U_ref); @@ -195,7 +197,7 @@ int main() //---------------------------------------------------------------------- size_t number_of_iterations = 0; int screen_output_interval = 100; - Real end_time = 1.0; + Real end_time = 5.0; Real output_interval = 0.02; /**< Time stamps for output of body states. */ Real dt = 0.0; /**< Default acoustic time step sizes. */ /** statistics for computing CPU time. */ @@ -230,6 +232,7 @@ int main() update_water_density_by_summation.exec(); update_air_density_by_summation.exec(); air_transport_correction.exec(); + water_transport_correction.exec(); air_viscous_acceleration.exec(); water_viscous_acceleration.exec(); diff --git a/tests/2d_examples/test_2d_taylor_green/taylor_green.cpp b/tests/2d_examples/test_2d_taylor_green/taylor_green.cpp index 0c54f4bcbf..8595555f97 100644 --- a/tests/2d_examples/test_2d_taylor_green/taylor_green.cpp +++ b/tests/2d_examples/test_2d_taylor_green/taylor_green.cpp @@ -100,7 +100,7 @@ int main(int ac, char *av[]) Dynamics1Level density_relaxation(water_block_inner); InteractionWithUpdate update_density_by_summation(water_block_inner); InteractionDynamics viscous_acceleration(water_block_inner); - InteractionDynamics transport_velocity_correction(water_block_inner); + InteractionDynamics> transport_velocity_correction(water_block_inner); ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); ReduceDynamics get_fluid_time_step_size(water_block); PeriodicConditionUsingCellLinkedList periodic_condition_x(water_block, water_block.getBodyShapeBounds(), xAxis); diff --git a/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp b/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp index 9a24f9b250..aa0b5956ca 100644 --- a/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp +++ b/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp @@ -315,7 +315,7 @@ int main(int ac, char *av[]) /** Computing viscous acceleration. */ InteractionDynamics viscous_acceleration(water_block_complex); /** Impose transport velocity formulation. */ - InteractionDynamics transport_velocity_correction(water_block_complex); + InteractionDynamics> transport_velocity_correction(water_block_complex); /** Computing vorticity in the flow. */ InteractionDynamics compute_vorticity(water_block_inner); /** Inflow boundary condition. */ diff --git a/tests/2d_examples/test_2d_throat/throat.cpp b/tests/2d_examples/test_2d_throat/throat.cpp index b17a389bee..bcbe226713 100644 --- a/tests/2d_examples/test_2d_throat/throat.cpp +++ b/tests/2d_examples/test_2d_throat/throat.cpp @@ -174,7 +174,7 @@ int main(int ac, char *av[]) InteractionSplit> implicit_viscous_damping(fluid_block_complex, "Velocity", mu_f); // impose transport velocity - InteractionDynamics transport_velocity_correction(fluid_block_complex); + InteractionDynamics> transport_velocity_correction(fluid_block_complex); // computing vorticity in the flow InteractionDynamics compute_vorticity(fluid_block_inner); //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_two_phase_dambreak/two_phase_dambreak.cpp b/tests/2d_examples/test_2d_two_phase_dambreak/two_phase_dambreak.cpp index 027f5b77a4..e6402e3402 100644 --- a/tests/2d_examples/test_2d_two_phase_dambreak/two_phase_dambreak.cpp +++ b/tests/2d_examples/test_2d_two_phase_dambreak/two_phase_dambreak.cpp @@ -60,7 +60,7 @@ int main() update_water_density_by_summation(water_wall_contact, water_air_complex.getInnerRelation()); InteractionWithUpdate update_air_density_by_summation(air_wall_contact, air_water_complex); - InteractionDynamics + InteractionDynamics> air_transport_correction(air_wall_contact, air_water_complex); /** Time step size without considering sound wave speed. */ ReduceDynamics get_water_advection_time_step_size(water_block, U_ref); diff --git a/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp b/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp index 4327321db1..0b6ac2244d 100644 --- a/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp +++ b/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp @@ -54,8 +54,10 @@ int main() InteractionWithUpdate update_air_density_by_summation(air_wall_contact, air_water_complex); /** transport formulation for regularizing particle distribution. */ - InteractionDynamics + InteractionDynamics> air_transport_correction(air_wall_contact, air_water_complex, 0.05); + InteractionDynamics> + water_transport_correction(water_wall_contact, water_air_complex, 0.05); /** Time step size without considering sound wave speed. */ ReduceDynamics get_water_advection_time_step_size(water_block, U_ref); ReduceDynamics get_air_advection_time_step_size(air_block, U_ref); @@ -135,6 +137,7 @@ int main() update_water_density_by_summation.exec(); update_air_density_by_summation.exec(); air_transport_correction.exec(); + water_transport_correction.exec(); air_viscous_acceleration.exec(); water_viscous_acceleration.exec(); diff --git a/tests/3d_examples/test_3d_stfb/stfb.cpp b/tests/3d_examples/test_3d_stfb/stfb.cpp index a61ec2684a..b9ebbe70f0 100644 --- a/tests/3d_examples/test_3d_stfb/stfb.cpp +++ b/tests/3d_examples/test_3d_stfb/stfb.cpp @@ -76,7 +76,7 @@ int main(int ac, char *av[]) InteractionWithUpdate free_stream_surface_indicator(water_block_complex); /** Impose transport velocity formulation. */ - InteractionDynamics + InteractionDynamics>> transport_velocity_correction(water_block_complex); /*-------------------------------------------------------------------------------*/ //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_stlw/stlw.cpp b/tests/3d_examples/test_3d_stlw/stlw.cpp index 9885fb0d01..12c548b98e 100644 --- a/tests/3d_examples/test_3d_stlw/stlw.cpp +++ b/tests/3d_examples/test_3d_stlw/stlw.cpp @@ -56,7 +56,7 @@ int main(int ac, char *av[]) InteractionWithUpdate free_stream_surface_indicator(water_block_complex); /** Impose transport velocity formulation. */ - InteractionDynamics + InteractionDynamics>> transport_velocity_correction(water_block_complex); /*-------------------------------------------------------------------------------*/ //---------------------------------------------------------------------- diff --git a/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.cpp b/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.cpp index 193642cd06..c0cebe7b5a 100644 --- a/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.cpp +++ b/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.cpp @@ -86,7 +86,7 @@ NonReflectiveBoundaryVariableCorrection::NonReflectiveBoundaryVariableCorrection fluid_(DynamicCast(this, particles_->getBaseMaterial())), rho_farfield_(0.0), sound_speed_(0.0), vel_farfield_(Vecd::Zero()), rho_(particles_->rho_), p_(*particles_->getVariableByName("Pressure")), Vol_(particles_->Vol_), vel_(particles_->vel_), mom_(*particles_->getVariableByName("Momentum")), pos_(particles_->pos_), - surface_indicator_(*particles_->getVariableByName("SurfaceIndicator")) + indicator_(*particles_->getVariableByName("Indicator")) { particles_->registerVariable(n_, "NormalDirection"); particles_->registerVariable(inner_weight_summation_, "InnerWeightSummation"); @@ -99,7 +99,7 @@ NonReflectiveBoundaryVariableCorrection::NonReflectiveBoundaryVariableCorrection //=================================================================================================// void NonReflectiveBoundaryVariableCorrection::initialization(size_t index_i, Real dt) { - if (surface_indicator_[index_i] == 1) + if (indicator_[index_i] == 1) { const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) @@ -114,7 +114,7 @@ void NonReflectiveBoundaryVariableCorrection::initialization(size_t index_i, Rea void NonReflectiveBoundaryVariableCorrection::interaction(size_t index_i, Real dt) { Shape &body_shape = *sph_body_.body_shape_; - if (surface_indicator_[index_i] == 1 || surface_inner_particle_indicator_[index_i] == 1) + if (indicator_[index_i] == 1 || surface_inner_particle_indicator_[index_i] == 1) { Vecd normal_direction = body_shape.findNormalDirection(pos_[index_i]); n_[index_i] = normal_direction; @@ -133,7 +133,7 @@ void NonReflectiveBoundaryVariableCorrection::interaction(size_t index_i, Real d for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - if (surface_indicator_[index_j] != 1) + if (indicator_[index_j] != 1) { Real W_ij = inner_neighborhood.W_ij_[n]; inner_weight_summation_[index_i] += W_ij * Vol_[index_j]; @@ -159,7 +159,7 @@ void NonReflectiveBoundaryVariableCorrection::interaction(size_t index_i, Real d for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - if (surface_indicator_[index_j] != 1) + if (indicator_[index_j] != 1) { rho_summation += rho_[index_j]; vel_summation += vel_[index_j]; @@ -182,7 +182,7 @@ void NonReflectiveBoundaryVariableCorrection::interaction(size_t index_i, Real d for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - if (surface_indicator_[index_j] != 1) + if (indicator_[index_j] != 1) { Real W_ij = inner_neighborhood.W_ij_[n]; inner_weight_summation_[index_i] += W_ij * Vol_[index_j]; @@ -203,7 +203,7 @@ void NonReflectiveBoundaryVariableCorrection::interaction(size_t index_i, Real d void NonReflectiveBoundaryVariableCorrection::update(size_t index_i, Real dt) { Shape &body_shape = *sph_body_.body_shape_; - if (surface_indicator_[index_i] == 1 || surface_inner_particle_indicator_[index_i] == 1) + if (indicator_[index_i] == 1 || surface_inner_particle_indicator_[index_i] == 1) { Vecd normal_direction = body_shape.findNormalDirection(pos_[index_i]); n_[index_i] = normal_direction; diff --git a/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.h b/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.h index e4687c59b5..e1468a512f 100644 --- a/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.h +++ b/tests/user_examples/extra_src/shared/common_weakly_compressible_eulerian_classes.h @@ -281,7 +281,7 @@ class NonReflectiveBoundaryVariableCorrection : public LocalDynamics, public Dat StdLargeVec n_; StdLargeVec inner_weight_summation_, rho_average_, vel_normal_average_; StdLargeVec vel_tangential_average_, vel_average_; - StdLargeVec &surface_indicator_; + StdLargeVec &indicator_; StdLargeVec surface_inner_particle_indicator_; }; } // namespace SPH diff --git a/tests/user_examples/test_2d_eulerian_flow_around_cylinder_LG/2d_eulerian_flow_around_cylinder_LG.cpp b/tests/user_examples/test_2d_eulerian_flow_around_cylinder_LG/2d_eulerian_flow_around_cylinder_LG.cpp index a344bb604e..9a85e581bf 100644 --- a/tests/user_examples/test_2d_eulerian_flow_around_cylinder_LG/2d_eulerian_flow_around_cylinder_LG.cpp +++ b/tests/user_examples/test_2d_eulerian_flow_around_cylinder_LG/2d_eulerian_flow_around_cylinder_LG.cpp @@ -34,7 +34,7 @@ int main(int ac, char *av[]) (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) ? water_block.generateParticles(io_environment, water_block.getName()) : water_block.generateParticles(); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); SolidBody cylinder(sph_system, makeShared("Cylinder")); cylinder.defineAdaptationRatios(1.3, 2.0); diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp index d03dc95bbc..353faf59ad 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp +++ b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp @@ -113,7 +113,7 @@ int main(int ac, char *av[]) InteractionWithUpdate update_fluid_density(water_block_complex); /** We can output a method-specific particle data for debug */ water_block.addBodyStateForRecording("Pressure"); - water_block.addBodyStateForRecording("SurfaceIndicator"); + water_block.addBodyStateForRecording("Indicator"); /** Time step size without considering sound wave speed. */ ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); /** Time step size with considering sound wave speed. */ @@ -128,7 +128,7 @@ int main(int ac, char *av[]) /** Computing viscous acceleration. */ InteractionDynamics viscous_acceleration(water_block_complex); /** Impose transport velocity formulation. */ - InteractionDynamics transport_velocity_correction(water_block_complex); + InteractionDynamics>> transport_velocity_correction(water_block_complex); /** Computing vorticity in the flow. */ InteractionDynamics compute_vorticity(water_block_inner); //----------------------------------------------------------------------