From 9f48fb1af66d228031caa012b5871d5ceb27bfa7 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 00:39:23 +0200 Subject: [PATCH 1/7] add plastic materials (soils) --- .../shared/general_continuum/all_continuum.h | 1 + .../continuum_dynamics_complex.h | 58 ++++ .../continuum_dynamics_complex.hpp | 135 ++++++++ .../continuum_dynamics_inner.cpp | 71 ++++- .../continuum_dynamics_inner.h | 99 +++++- .../continuum_dynamics_inner.hpp | 139 +++++++++ .../general_continuum/continuum_particles.cpp | 48 +++ .../general_continuum/continuum_particles.h | 30 ++ .../general_continuum/general_continuum.cpp | 62 ++++ .../general_continuum/general_continuum.h | 31 ++ .../riemann_solver_extra.cpp | 14 + .../general_continuum/riemann_solver_extra.h | 26 ++ ...etting_coupled_spatial_temporal_method.cpp | 29 -- .../wetting_coupled_spatial_temporal_method.h | 107 ------- .../test_2d_column_collapse/CMakeLists.txt | 21 ++ .../column_collapse.cpp | 278 +++++++++++++++++ .../GranularBody_TotalMechanicalEnergy.xml | 36 +++ ...ody_TotalMechanicalEnergy_Run_0_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_1_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_2_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_3_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_4_result.xml | 9 + ...Body_TotalMechanicalEnergy_dtwdistance.xml | 4 + ...larBody_TotalMechanicalEnergy_runtimes.dat | 3 + .../regression_test_tool.py | 37 +++ .../test_3d_repose_angle/CMakeLists.txt | 21 ++ .../FluidObserver_Pressure_runtimes.dat | 3 + .../GranularBody_TotalMechanicalEnergy.xml | 21 ++ ...ody_TotalMechanicalEnergy_Run_0_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_1_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_2_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_3_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_4_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_5_result.xml | 9 + ...ody_TotalMechanicalEnergy_Run_6_result.xml | 9 + ...Body_TotalMechanicalEnergy_dtwdistance.xml | 4 + ...larBody_TotalMechanicalEnergy_runtimes.dat | 3 + .../regression_test_tool.py | 38 +++ .../test_3d_repose_angle/repose_angle.cpp | 295 ++++++++++++++++++ 39 files changed, 1583 insertions(+), 139 deletions(-) create mode 100644 tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp create mode 100644 tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.h delete mode 100644 tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp delete mode 100644 tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h create mode 100644 tests/user_examples/test_2d_column_collapse/CMakeLists.txt create mode 100644 tests/user_examples/test_2d_column_collapse/column_collapse.cpp create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat create mode 100644 tests/user_examples/test_2d_column_collapse/regression_test_tool/regression_test_tool.py create mode 100644 tests/user_examples/test_3d_repose_angle/CMakeLists.txt create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/FluidObserver_Pressure_runtimes.dat create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat create mode 100644 tests/user_examples/test_3d_repose_angle/regression_test_tool/regression_test_tool.py create mode 100644 tests/user_examples/test_3d_repose_angle/repose_angle.cpp diff --git a/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h b/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h index 258ea393d6..04bd078b74 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h +++ b/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h @@ -34,3 +34,4 @@ #include "all_continuum_dynamics.h" #include "general_continuum.h" #include "continuum_particles.h" +#include "riemann_solver_extra.h" diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h index 6db8a9da24..c0e82470ab 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h @@ -17,6 +17,11 @@ typedef DataDelegateContact FSIContactData; +typedef DataDelegateContact +PlasticContinuumWallData; +typedef DataDelegateContact +PlasticContinuumContactData; +typedef DataDelegateContact PlasticFSIContactData; /** * @class BaseShearStressRelaxation1stHalfType */ @@ -51,6 +56,59 @@ class BaseShearStressRelaxation2ndHalfWithWall : public fluid_dynamics::Interact }; using ShearStressRelaxation2ndHalfWithWall = BaseShearStressRelaxation2ndHalfWithWall; + +//=============================================Plasticity==========================================// +//=================================================================================================// +//==============================BaseStressRelaxation1stHalfWithWall=================================// +//=================================================================================================// +template +class BaseStressRelaxation1stHalfWithWall : public fluid_dynamics::InteractionWithWall +{ +public: + template + BaseStressRelaxation1stHalfWithWall(Args &&...args) + : fluid_dynamics::InteractionWithWall(std::forward(args)...) {}; + virtual ~BaseStressRelaxation1stHalfWithWall() {}; + void interaction(size_t index_i, Real dt = 0.0); + +protected: + virtual Vecd computeNonConservativeAcceleration(size_t index_i) override; +}; +using StressRelaxation1stHalfWithWall = BaseStressRelaxation1stHalfWithWall; +using StressRelaxation1stHalfRiemannWithWall = BaseStressRelaxation1stHalfWithWall; +using StressRelaxation1stHalfDissipativeRiemannfWithWall = BaseStressRelaxation1stHalfWithWall; +//=================================================================================================// +//===============================BaseStressRelaxation2ndHalfWithWall======================================// +//=================================================================================================// + +template +class BaseStressRelaxation2ndHalfWithWall : public fluid_dynamics::InteractionWithWall +{ +public: + template + BaseStressRelaxation2ndHalfWithWall(Args &&...args) + : fluid_dynamics::InteractionWithWall(std::forward(args)...) {}; + virtual ~BaseStressRelaxation2ndHalfWithWall() {}; + void interaction(size_t index_i, Real dt = 0.0); +}; +using StressRelaxation2ndHalfWithWall = BaseStressRelaxation2ndHalfWithWall; +using StressRelaxation2ndHalfRiemannWithWall = BaseStressRelaxation2ndHalfWithWall; +using StressRelaxation2ndHalfDissipativeRiemannWithWall = BaseStressRelaxation2ndHalfWithWall; + +/** +* @class StressDiffusionWithWall +*/ +template +class BaseStressDiffusionWithWall : public fluid_dynamics::InteractionWithWall +{ +public: + template + BaseStressDiffusionWithWall(Args &&...args) + : fluid_dynamics::InteractionWithWall(std::forward(args)...) {}; + virtual ~BaseStressDiffusionWithWall() {}; + void interaction(size_t index_i, Real dt = 0.0); +}; +using StressDiffusionWithWall = BaseStressDiffusionWithWall; } // namespace continuum_dynamics } // namespace SPH diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp index 37e7286f35..6f86f35f56 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp @@ -70,5 +70,140 @@ void BaseShearStressRelaxation2ndHalfWithWallvelocity_gradient_[index_i] += velocity_gradient; } + +//============================================Plasticity===========================================// +//=================================================================================================// +//============BaseStressRelaxation1stHalfWithWall & BaseStressRelaxation2ndHalfWithWall============// +//=================================================================================================// +template +Vecd BaseStressRelaxation1stHalfWithWall::computeNonConservativeAcceleration(size_t index_i) +{ + return this->acc_prior_[index_i]; +} + +template +void BaseStressRelaxation1stHalfWithWall::interaction(size_t index_i, Real dt) +{ + BaseStressRelaxation1stHalfType::interaction(index_i, dt); + + Vecd acc_prior_i = computeNonConservativeAcceleration(index_i); + + //Vecd acceleration = Vecd::Zero(); + Vecd acceleration = acc_prior_i; + Real rho_dissipation(0); + + Matd stress_tensor_i = this->reduceTensor(this->stress_tensor_3D_[index_i]); + //Real rho_i = this->rho_[index_i]; + //Real rho_in_wall = this->plastic_continuum_.getDensity(); + + for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) + { + StdLargeVec& acc_ave_k = *(this->wall_acc_ave_[k]); + Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != wall_neighborhood.current_size_; ++n) + { + size_t index_j = wall_neighborhood.j_[n]; + Vecd& e_ij = wall_neighborhood.e_ij_[n]; + Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n]; + Real r_ij = wall_neighborhood.r_ij_[n]; + Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; + + Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); + Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_wall_external_acceleration); + //Real p_in_wall = this->p_[index_i]; + Matd stress_tensor_in_wall = stress_tensor_i; + //acceleration += rho_in_wall * (stress_tensor_i / (rho_i * rho_i) + stress_tensor_in_wall / (rho_in_wall * rho_in_wall)) * nablaW_ijV_j; + acceleration += (stress_tensor_i + stress_tensor_in_wall) * nablaW_ijV_j; + + rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; + } + } + //this->acc_[index_i] += acceleration; + this->acc_[index_i] += acceleration / this->rho_[index_i]; + this->drho_dt_[index_i] += rho_dissipation * this->rho_[index_i]; +} +//=================================================================================================// +template +void BaseStressRelaxation2ndHalfWithWall::interaction(size_t index_i, Real dt) +{ + BaseStressRelaxation2ndHalfType::interaction(index_i, dt); + + Real density_change_rate = 0.0; + Vecd p_dissipation = Vecd::Zero(); + + Vecd vel_i = this->vel_[index_i]; + Matd velocity_gradient = Matd::Zero(); + + for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) + { + StdLargeVec& vel_ave_k = *(this->wall_vel_ave_[k]); + StdLargeVec& n_k = *(this->wall_n_[k]); + + Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != wall_neighborhood.current_size_; ++n) + { + size_t index_j = wall_neighborhood.j_[n]; + Vecd& e_ij = wall_neighborhood.e_ij_[n]; + Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n]; + Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; + + Real beta = 1.9; + Vecd vel_in_wall = (1 - beta) * vel_i + beta * vel_ave_k[index_j]; + + Matd velocity_gradient_ij = -(vel_i - vel_in_wall) * nablaW_ijV_j.transpose(); + velocity_gradient += velocity_gradient_ij; + + density_change_rate += (this->vel_[index_i] - vel_in_wall).dot(e_ij) * dW_ijV_j; + // Real u_jump = 2.0 * (this->vel_[index_i] - vel_ave_k[index_j]).dot(n_k[index_j]); + // p_dissipation += this->riemann_solver_.DissipativePJump(u_jump) * dW_ijV_j * n_k[index_j]; + Vecd u_jump = 2.0 * (this->vel_[index_i] - vel_ave_k[index_j]); + p_dissipation += this->riemann_solver_.DissipativePJumpExtra(u_jump, n_k[index_j]) * dW_ijV_j; + // Real temp = u_jump.dot(n_k[index_j]); + // p_dissipation += this->riemann_solver_.DissipativePJump(temp) * dW_ijV_j * e_ij; + } + } + this->drho_dt_[index_i] += density_change_rate * this->rho_[index_i]; + this->velocity_gradient_[index_i] += velocity_gradient; + this->acc_[index_i] += p_dissipation / this->rho_[index_i]; +} + +//=================================================================================================// + //================================BaseStressDiffusionWithWall======================================// + //=================================================================================================// +template +void BaseStressDiffusionWithWall::interaction(size_t index_i, Real dt) +{ + BaseStressDiffusionType::interaction(index_i, dt); + + Real rho_in_wall = this->plastic_continuum_.getDensity(); + Real gravity = abs(this->acc_prior_[index_i](1, 0)); + Real density = rho_in_wall; + Mat3d diffusion_stress_rate_ = Mat3d::Zero(); + Mat3d diffusion_stress_ = Mat3d::Zero(); + Mat3d stress_tensor_i = this->stress_tensor_3D_[index_i]; + for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) + { + Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != wall_neighborhood.current_size_; ++n) + { + size_t index_j = wall_neighborhood.j_[n]; + Real r_ij = wall_neighborhood.r_ij_[n]; + Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n]; + + Real y_ij = this->pos_[index_i](1, 0) - this->pos_[index_j](1, 0); + //stress boundary condition + Mat3d stress_tensor_j = stress_tensor_i; + + diffusion_stress_ = stress_tensor_i - stress_tensor_j; + diffusion_stress_(0, 0) = diffusion_stress_(0, 0) - (1 - sin(this->fai_)) * density * gravity * y_ij; + diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij; + diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(this->fai_)) * density * gravity * y_ij; + + diffusion_stress_rate_ += 2 * this->zeta_ * this->smoothing_length_ * this->sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * this->smoothing_length_); + } + } + this->stress_rate_3D_[index_i] += diffusion_stress_rate_; +} + } // namespace continuum_dynamics } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp index cc91c4bdcd..39a96314ed 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp @@ -5,8 +5,8 @@ namespace SPH namespace continuum_dynamics { ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body) - : LocalDynamics(sph_body), ContinuumDataSimple(sph_body), - pos_(particles_->pos_), vel_(particles_->vel_) {} + : LocalDynamics(sph_body), PlasticContinuumDataSimple(sph_body), + pos_(particles_->pos_), vel_(particles_->vel_), stress_tensor_3D_(particles_->stress_tensor_3D_) {} //=================================================================================================// ContinuumAcousticTimeStepSize::ContinuumAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL) : fluid_dynamics::AcousticTimeStepSize(sph_body, acousticCFL) {} @@ -282,5 +282,72 @@ namespace SPH { vel_[index_i] -= velocity_correction_; } + //=================================================================================================// + //================================================Plastic==========================================// + //=================================================================================================// + BaseRelaxationPlastic::BaseRelaxationPlastic(BaseInnerRelation& inner_relation) + : LocalDynamics(inner_relation.getSPHBody()), PlasticContinuumDataInner(inner_relation), + plastic_continuum_(particles_->plastic_continuum_), rho_(particles_->rho_), + p_(*particles_->getVariableByName("Pressure")), drho_dt_(*particles_->registerSharedVariable("DensityChangeRate")), pos_(particles_->pos_), vel_(particles_->vel_), acc_(particles_->acc_), acc_prior_(particles_->acc_prior_), + stress_tensor_3D_(particles_->stress_tensor_3D_), strain_tensor_3D_(particles_->strain_tensor_3D_), + stress_rate_3D_(particles_->stress_rate_3D_), strain_rate_3D_(particles_->strain_rate_3D_), + elastic_strain_tensor_3D_(particles_->elastic_strain_tensor_3D_), elastic_strain_rate_3D_(particles_->elastic_strain_rate_3D_) {} + Matd BaseRelaxationPlastic::reduceTensor(Mat3d tensor_3d) + { + Matd tensor_2d; + for (int i = 0; i < (Real)Dimensions; i++) + { + for (int j = 0; j < (Real)Dimensions; j++) + { + tensor_2d(i, j) = tensor_3d(i, j); + } + } + return tensor_2d; + } + + Mat3d BaseRelaxationPlastic::increaseTensor(Matd tensor_2d) + { + Mat3d tensor_3d = Mat3d::Zero(); + for (int i = 0; i < (Real)Dimensions; i++) + { + for (int j = 0; j < (Real)Dimensions; j++) + { + tensor_3d(i, j) = tensor_2d(i, j); + } + } + return tensor_3d; + } + + //====================================================================================// + //===============================StressDiffusion======================================// + //====================================================================================// + StressDiffusion::StressDiffusion(BaseInnerRelation& inner_relation) + : BaseRelaxationPlastic(inner_relation), fai_(DynamicCast(this, plastic_continuum_).getFrictionAngle()), smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), + sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {} + void StressDiffusion::interaction(size_t index_i, Real dt) + { + Real gravity = abs(acc_prior_[index_i](1, 0)); + Real density = plastic_continuum_.getDensity(); + Mat3d diffusion_stress_rate_ = Mat3d::Zero(); + Mat3d diffusion_stress_ = Mat3d::Zero(); + + 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]; + Real r_ij = inner_neighborhood.r_ij_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; + Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); + //calculate Psi, equation (34) in Feng_2021_CG + diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; + diffusion_stress_(0, 0) = diffusion_stress_(0, 0) - (1 - sin(fai_)) * density * gravity * y_ij; + diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij; + diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(fai_)) * density * gravity * y_ij; + + //diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_ * smoothing_length_); + diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); + } + stress_rate_3D_[index_i] = diffusion_stress_rate_; + } } // namespace continuum_dynamics } // namespace SPH diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h index 9f9919a62f..2ccc7d8714 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h @@ -4,6 +4,7 @@ #include "fluid_dynamics_inner.hpp" #include "general_continuum.h" #include "continuum_particles.h" +#include "riemann_solver_extra.h" namespace SPH { @@ -12,7 +13,10 @@ namespace continuum_dynamics typedef DataDelegateSimple ContinuumDataSimple; typedef DataDelegateInner ContinuumDataInner; -class ContinuumInitialCondition : public LocalDynamics, public ContinuumDataSimple +typedef DataDelegateSimple PlasticContinuumDataSimple; +typedef DataDelegateInner PlasticContinuumDataInner; + +class ContinuumInitialCondition : public LocalDynamics, public PlasticContinuumDataSimple { public: explicit ContinuumInitialCondition(SPHBody &sph_body); @@ -20,6 +24,7 @@ class ContinuumInitialCondition : public LocalDynamics, public ContinuumDataSimp protected: StdLargeVec &pos_, &vel_; + StdLargeVec& stress_tensor_3D_; }; class ContinuumAcousticTimeStepSize : public fluid_dynamics::AcousticTimeStepSize @@ -254,6 +259,98 @@ class ConstrainSolidBodyMassCenter : public LocalDynamics, public ContinuumDataS void update(size_t index_i, Real dt = 0.0); }; + +//=================================================================================================// +//================================================Plastic==========================================// +//=================================================================================================// + +class BaseRelaxationPlastic : public LocalDynamics, public PlasticContinuumDataInner +{ +public: + explicit BaseRelaxationPlastic(BaseInnerRelation& inner_relation); + virtual ~BaseRelaxationPlastic() {}; + Matd reduceTensor(Mat3d tensor_3d); + Mat3d increaseTensor(Matd tensor_2d); + +protected: + PlasticContinuum& plastic_continuum_; + StdLargeVec& rho_, & p_, & drho_dt_; + StdLargeVec& pos_, & vel_, & acc_, & acc_prior_; + StdLargeVec& stress_tensor_3D_, & strain_tensor_3D_, & stress_rate_3D_, & strain_rate_3D_; + StdLargeVec& elastic_strain_tensor_3D_, & elastic_strain_rate_3D_; + +}; + +//=================================================================================================// + //===============================BaseStressRelaxation1stHalf======================================// + //=================================================================================================// + /** + * @class BaseIntegration1stHalf + * @brief Template class for pressure relaxation scheme with the Riemann solver + * as template variable + */ +template +class BaseStressRelaxation1stHalf : public BaseRelaxationPlastic +{ +public: + explicit BaseStressRelaxation1stHalf(BaseInnerRelation& inner_relation); + virtual ~BaseStressRelaxation1stHalf() {}; + RiemannSolverType riemann_solver_; + void initialization(size_t index_i, Real dt = 0.0); + void interaction(size_t index_i, Real dt = 0.0); + void update(size_t index_i, Real dt = 0.0); + +protected: + virtual Vecd computeNonConservativeAcceleration(size_t index_i); + + StdLargeVec& velocity_gradient_; +}; +using StressRelaxation1stHalf = BaseStressRelaxation1stHalf; +using StressRelaxation1stHalfRiemann = BaseStressRelaxation1stHalf; +using StressRelaxation1stHalfDissipativeRiemann = BaseStressRelaxation1stHalf; + +//=================================================================================================// +//===============================BaseStressRelaxation2ndHalf======================================// +//=================================================================================================// +/** + * @class BaseIntegration2ndHalf + * @brief Template density relaxation scheme with different Riemann solver + */ +template +class BaseStressRelaxation2ndHalf : public BaseRelaxationPlastic +{ +public: + explicit BaseStressRelaxation2ndHalf(BaseInnerRelation& inner_relation); + virtual ~BaseStressRelaxation2ndHalf() {}; + RiemannSolverType riemann_solver_; + void initialization(size_t index_i, Real dt = 0.0); + void interaction(size_t index_i, Real dt = 0.0); + void update(size_t index_i, Real dt = 0.0); + +protected: + + StdLargeVec& velocity_gradient_; + StdLargeVec& Vol_, & mass_, & von_mises_stress_; + StdLargeVec& acc_deviatoric_plastic_strain_, & vertical_stress_; + Real E_, nu_; + +}; +using StressRelaxation2ndHalf = BaseStressRelaxation2ndHalf; +using StressRelaxation2ndHalfRiemann = BaseStressRelaxation2ndHalf; +using StressRelaxation2ndHalfDissipativeRiemann = BaseStressRelaxation2ndHalf; +/** +* @class StressDiffusion +*/ +class StressDiffusion : public BaseRelaxationPlastic +{ +public: + explicit StressDiffusion(BaseInnerRelation& inner_relation); + virtual ~StressDiffusion() {}; + void interaction(size_t index_i, Real dt = 0.0); +protected: + Real zeta_ = 0.1, fai_; //diffusion coeficient + Real smoothing_length_, sound_speed_; +}; } // namespace continuum_dynamics } // namespace SPH #endif // CONTINUUM_DYNAMICS_INNER_H \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp index a653c918c6..cb626b7f93 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp @@ -24,5 +24,144 @@ namespace SPH { this->vel_[index_i] += (this->acc_prior_[index_i] + this->acc_[index_i] + this->acc_shear_[index_i]) * dt; } + + //=================================================================================================// + //================Plastic: BaseStressRelaxation1stHalf & BaseStressRelaxation2ndHalf==============// + //=================================================================================================// + template + BaseStressRelaxation1stHalf::BaseStressRelaxation1stHalf(BaseInnerRelation& inner_relation) + : BaseRelaxationPlastic(inner_relation), riemann_solver_(plastic_continuum_, plastic_continuum_), + velocity_gradient_(particles_->velocity_gradient_) {} + //=================================================================================================// + template + void BaseStressRelaxation1stHalf::initialization(size_t index_i, Real dt) + { + rho_[index_i] += drho_dt_[index_i] * dt * 0.5; + p_[index_i] = plastic_continuum_.getPressure(rho_[index_i]); + //pressure_replace_by_volumn_part + //p_[index_i] = - stress_tensor_3D_[index_i].trace() /3; + + pos_[index_i] += vel_[index_i] * dt * 0.5; + + stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt * 0.5; + strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt * 0.5; + + //calculate elastic strain + elastic_strain_tensor_3D_[index_i] += elastic_strain_rate_3D_[index_i] * dt * 0.5; + } + //=================================================================================================// + template + void BaseStressRelaxation1stHalf::update(size_t index_i, Real dt) + { + vel_[index_i] += (acc_prior_[index_i] + acc_[index_i]) * dt; + } + //=================================================================================================// + template + Vecd BaseStressRelaxation1stHalf::computeNonConservativeAcceleration(size_t index_i) + { + Vecd acceleration = acc_prior_[index_i] * rho_[index_i]; + const 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]; + Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; + const Vecd& e_ij = inner_neighborhood.e_ij_[n]; + + acceleration += (p_[index_i] - p_[index_j]) * dW_ijV_j * e_ij; + } + return acceleration / rho_[index_i]; + } + //=================================================================================================// + template + void BaseStressRelaxation1stHalf::interaction(size_t index_i, Real dt) + { + Vecd acceleration = Vecd::Zero(); + Real rho_dissipation(0); + Real rho_i = rho_[index_i]; + Matd stress_tensor_i = reduceTensor(stress_tensor_3D_[index_i]); + const 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]; + Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; + Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + Matd stress_tensor_j = reduceTensor(stress_tensor_3D_[index_j]); + //acceleration += rho_[index_j] * (stress_tensor_i / (rho_i * rho_i) + stress_tensor_j / (rho_[index_j] * rho_[index_j])) * nablaW_ijV_j; + acceleration += rho_[index_j] * ((stress_tensor_i + stress_tensor_j) / (rho_i * rho_[index_j])) * nablaW_ijV_j; + rho_dissipation += riemann_solver_.DissipativeUJump(p_[index_i] - p_[index_j]) * dW_ijV_j; + } + acc_[index_i] += acceleration; + drho_dt_[index_i] = rho_dissipation * rho_[index_i]; + } + //=================================================================================================// + template + BaseStressRelaxation2ndHalf::BaseStressRelaxation2ndHalf(BaseInnerRelation& inner_relation) + : BaseRelaxationPlastic(inner_relation), riemann_solver_(plastic_continuum_, plastic_continuum_), + velocity_gradient_(particles_->velocity_gradient_), Vol_(particles_->Vol_), mass_(particles_->mass_), von_mises_stress_(particles_->von_mises_stress_), acc_deviatoric_plastic_strain_(particles_->acc_deviatoric_plastic_strain_), vertical_stress_(particles_->vertical_stress_), E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) {} + //=================================================================================================// + template + void BaseStressRelaxation2ndHalf::initialization(size_t index_i, Real dt) + { + pos_[index_i] += vel_[index_i] * dt * 0.5; + } + //=================================================================================================// + template + void BaseStressRelaxation2ndHalf::interaction(size_t index_i, Real dt) + { + Real density_change_rate(0); + Vecd p_dissipation = Vecd::Zero(); + Matd velocity_gradient = Matd::Zero(); + const 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]; + const Vecd& e_ij = inner_neighborhood.e_ij_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; + Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + + Vecd u_jump = vel_[index_i] - vel_[index_j]; + density_change_rate += u_jump.dot(e_ij) * dW_ijV_j; + p_dissipation += riemann_solver_.DissipativePJumpExtra(u_jump, e_ij) * dW_ijV_j; + // Real temp = u_jump.dot(e_ij); + // p_dissipation += riemann_solver_.DissipativePJump(temp) * dW_ijV_j * e_ij; + + Matd velocity_gradient_ij = -(vel_[index_i] - vel_[index_j]) * nablaW_ijV_j.transpose(); + velocity_gradient += velocity_gradient_ij; + } + drho_dt_[index_i] += density_change_rate * rho_[index_i]; + acc_[index_i] = p_dissipation / rho_[index_i]; + velocity_gradient_[index_i] = velocity_gradient; + } + //=================================================================================================// + template + void BaseStressRelaxation2ndHalf::update(size_t index_i, Real dt) + { + rho_[index_i] += drho_dt_[index_i] * dt * 0.5; + Vol_[index_i] = mass_[index_i] / rho_[index_i]; + + Mat3d velocity_gradient = increaseTensor(velocity_gradient_[index_i]); + Mat3d stress_tensor_rate_3D_ = plastic_continuum_.ConstitutiveRelationZ(velocity_gradient, stress_tensor_3D_[index_i]); + //consider diffusion stress with += + //stress_rate_3D_[index_i] = stress_tensor_rate_3D_; + stress_rate_3D_[index_i] += stress_tensor_rate_3D_; + //update stress tensor + stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt * 0.5; + + //For plasticity + Mat3d stress_tensor_ = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); + stress_tensor_3D_[index_i] = stress_tensor_; + //added + vertical_stress_[index_i] = stress_tensor_3D_[index_i](1, 1); + //Accumulated deviatoric plastic strains + strain_rate_3D_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); + strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt * 0.5; + //calculate elastic strain + Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1 / 3) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity(); + Real hydrostatic_pressure = (1 / 3) * stress_tensor_3D_[index_i].trace(); + elastic_strain_tensor_3D_[index_i] = deviatoric_stress / (2 * plastic_continuum_.getShearModulus(E_, nu_)) + hydrostatic_pressure * Mat3d::Identity() / (9 * plastic_continuum_.getBulkModulus(E_, nu_)); + Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D_[index_i]; + acc_deviatoric_plastic_strain_[index_i] = particles_->getDeviatoricPlasticStrain(plastic_strain_tensor_3D); + } } } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp index c195c1f653..4c169dd08c 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp @@ -43,4 +43,52 @@ namespace SPH { return n_[i]; }); } + //=================================================================================================// + PlasticContinuumParticles:: + PlasticContinuumParticles(SPHBody& sph_body, PlasticContinuum* plastic_continuum) + : ContinuumParticles(sph_body, plastic_continuum), plastic_continuum_(*plastic_continuum) {} + //=================================================================================================// + void PlasticContinuumParticles::initializeOtherVariables() + { + ContinuumParticles::initializeOtherVariables(); + + registerVariable(elastic_strain_tensor_3D_, "ElasticStrainTensor3D"); + registerVariable(elastic_strain_rate_3D_, "ElasticStrainRate3D"); + registerVariable(strain_tensor_3D_, "StrainTensor3D"); + registerVariable(stress_tensor_3D_, "StressTensor3D"); + registerVariable(strain_rate_3D_, "StrainRate3D"); + registerVariable(stress_rate_3D_, "StressRate3D"); + registerVariable(shear_stress_3D_, "ShearStress3D"); + registerVariable(shear_strain_3D_, "ShearStrain3D"); + registerVariable(shear_stress_rate_3D_, "ShearStressRate3D"); + registerVariable(shear_strain_rate_3D_, "ShearStrainRate3D"); + registerVariable(vertical_stress_, "VerticalStress"); + registerVariable(acc_deviatoric_plastic_strain_, "AccDeviatoricPlasticStrain"); + //---------------------------------------------------------------------- + // register sortable particle data + //---------------------------------------------------------------------- + registerSortableVariable("ElasticStrainTensor3D"); + registerSortableVariable("ElasticStrainRate3D"); + registerSortableVariable("StrainTensor3D"); + registerSortableVariable("StressTensor3D"); + registerSortableVariable("StrainRate3D"); + registerSortableVariable("StressRate3D"); + + registerSortableVariable("ShearStress3D"); + registerSortableVariable("ShearStrain3D"); + registerSortableVariable("ShearStressRate3D"); + registerSortableVariable("ShearStrainRate3D"); + + registerSortableVariable("VerticalStress"); + registerSortableVariable("AccDeviatoricPlasticStrain"); + + } + //=================================================================================================// + Real PlasticContinuumParticles::getDeviatoricPlasticStrain(Mat3d& strain_tensor) + { + Mat3d deviatoric_strain_tensor = strain_tensor - (1 / (Real)Dimensions) * strain_tensor.trace() * Mat3d::Identity(); + Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum(); + return sqrt(sum * 2 / 3); + } + } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h index 6f5ed620a8..baa22c422c 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h @@ -35,6 +35,36 @@ namespace SPH virtual void initializeOtherVariables() override; virtual ContinuumParticles *ThisObjectPtr() override { return this; }; }; + + class PlasticContinuumParticles : public ContinuumParticles + { + public: + StdLargeVec elastic_strain_tensor_3D_; + StdLargeVec elastic_strain_rate_3D_; + + StdLargeVec strain_tensor_3D_; + StdLargeVec stress_tensor_3D_; + StdLargeVec strain_rate_3D_; + StdLargeVec stress_rate_3D_; + + StdLargeVec shear_stress_3D_; + StdLargeVec shear_strain_3D_; + StdLargeVec shear_stress_rate_3D_; + StdLargeVec shear_strain_rate_3D_; + + StdLargeVec vertical_stress_; + StdLargeVec acc_deviatoric_plastic_strain_; + + Real getDeviatoricPlasticStrain(Mat3d& strain_tensor); + + PlasticContinuum& plastic_continuum_; + + PlasticContinuumParticles(SPHBody& sph_body, PlasticContinuum* plastic_continuum); + virtual ~PlasticContinuumParticles() {}; + + virtual void initializeOtherVariables() override; + virtual ContinuumParticles* ThisObjectPtr() override { return this; }; + }; } // namespace SPH #endif // CONTINUUM_PARTICLES_H \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp index 63f34a2aac..d247123aaa 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp @@ -26,4 +26,66 @@ namespace SPH Matd stress_rate = 2.0 * G_ * deviatoric_strain_rate + shear_stress * (spin_rate.transpose()) + spin_rate * shear_stress; return stress_rate; } + //=============================================================================================// + //==============================PlasticContinuum===============================================// + //=============================================================================================// + Real PlasticContinuum::getPressure(Real rho) + { + return p0_ * (rho / rho0_ - 1.0); + } + + Real PlasticContinuum::getDPConstantsA(Real friction_angle) + { + return tan(friction_angle) / sqrt(9 + 12 * tan(friction_angle) * tan(friction_angle)); + } + + Real PlasticContinuum::getDPConstantsK(Real cohesion, Real friction_angle) + { + return 3 * cohesion / sqrt(9 + 12 * tan(friction_angle) * tan(friction_angle)); + } + + Mat3d PlasticContinuum::ConstitutiveRelationZ(Mat3d& velocity_gradient, Mat3d& stress_tensor) + { + Real dim = 3; + Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); + Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose()); + Mat3d deviatoric_strain_rate = strain_rate - (1.0 / dim) * strain_rate.trace() * Mat3d::Identity(); + //consider the elastic part + Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor; + //consider the plastic part + Real stress_tensor_I1 = stress_tensor.trace(); // first invariant of stress + Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / dim) * stress_tensor.trace() * Mat3d::Identity(); + Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); + Real f = sqrt(stress_tensor_J2) + alpha_fai_ * stress_tensor_I1 - k_c_; + Real lambda_dot_ = 0; + Mat3d g = Mat3d::Zero(); + if ((f >= TinyReal) && (stress_tensor_J2 > TinyReal)) + { + Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum(); + lambda_dot_ = (3 * alpha_fai_ * K_ * strain_rate.trace() + (G_ / sqrt(stress_tensor_J2)) * deviatoric_stress_times_strain_rate) / (27 * alpha_fai_ * K_ * sin(psi_) + G_); + g = lambda_dot_ * (9 * K_ * sin(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2))); + } + Mat3d stress_rate_temp = stress_rate_elastic - g; + return stress_rate_temp; + } + + Mat3d PlasticContinuum::ReturnMapping(Mat3d& stress_tensor) + { + Real dim = 3; + + Real stress_tensor_I1 = stress_tensor.trace(); + if (-alpha_fai_ * stress_tensor_I1 + k_c_ < 0) + { + stress_tensor -= (1.0 / dim) * (stress_tensor_I1 - k_c_ / alpha_fai_) * Mat3d::Identity(); + } + stress_tensor_I1 = stress_tensor.trace(); + Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / dim) * stress_tensor.trace() * Mat3d::Identity(); + Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); + if (-alpha_fai_ * stress_tensor_I1 + k_c_ < sqrt(stress_tensor_J2)) + { + Real r = (-alpha_fai_ * stress_tensor_I1 + k_c_) / (sqrt(stress_tensor_J2) + TinyReal); + stress_tensor = r * deviatoric_stress_tensor + (1.0 / dim) * stress_tensor_I1 * Mat3d::Identity(); + } + return stress_tensor; + } } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h index 0e016ba79f..853c4d2442 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h +++ b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h @@ -41,5 +41,36 @@ class GeneralContinuum : public WeaklyCompressibleFluid virtual GeneralContinuum *ThisObjectPtr() override { return this; }; }; + +class PlasticContinuum : public GeneralContinuum +{ +protected: + Real c_; /*< cohesion */ + Real fai_; /*< friction angle */ + Real psi_; /*< dilatancy angle */ + Real alpha_fai_; /*< Drucker¨CPragerˇŻs constants */ + Real k_c_; /*< Drucker¨CPragerˇŻs constants */ +public: + explicit PlasticContinuum(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio, Real friction_angle, Real cohesion=0, Real dilatancy=0) + : GeneralContinuum(rho0, c0, youngs_modulus, poisson_ratio), + c_(cohesion), fai_(friction_angle), psi_(dilatancy), alpha_fai_(0.0), k_c_(0.0) + { + material_type_name_ = "PlasticContinuum"; + alpha_fai_ = getDPConstantsA(friction_angle); + k_c_ = getDPConstantsK(cohesion, friction_angle); + }; + virtual ~PlasticContinuum() {}; + + + Real getDPConstantsA(Real friction_angle); + Real getDPConstantsK(Real cohesion, Real friction_angle); + Real getFrictionAngle() { return fai_; }; + virtual Real getPressure(Real rho); + + virtual Mat3d ConstitutiveRelationZ(Mat3d& velocity_gradient, Mat3d& stress_tensor); + virtual Mat3d ReturnMapping(Mat3d& stress_tensor); + + virtual GeneralContinuum* ThisObjectPtr() override { return this; }; +}; } // namespace SPH #endif // GENERAL_CONTINUUM_H diff --git a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp new file mode 100644 index 0000000000..06f6f6cba7 --- /dev/null +++ b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp @@ -0,0 +1,14 @@ +#include "riemann_solver_extra.h" + +namespace SPH +{ + Vecd AcousticRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) + { + return rho0c0_geo_ave_ * u_jump * SMIN(3.0 * SMAX(u_jump.dot(e_ij) * inv_c_ave_, 0.0), 1.0); + } + + Vecd DissipativeRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) + { + return rho0c0_geo_ave_ * u_jump; + } +} // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.h b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.h new file mode 100644 index 0000000000..c7c2b1b6bc --- /dev/null +++ b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.h @@ -0,0 +1,26 @@ +#ifndef RIEMANN_SOLVER_EXTRA_H +#define RIEMANN_SOLVER_EXTRA_H + +#include "riemann_solver.h" + +namespace SPH +{ + class AcousticRiemannSolverExtra : public AcousticRiemannSolver + { + public: + template + AcousticRiemannSolverExtra(FluidI& fluid_i, FluidJ& fluid_j) + : AcousticRiemannSolver(fluid_i, fluid_j) {}; + Vecd DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij); + }; + + class DissipativeRiemannSolverExtra : public DissipativeRiemannSolver + { + public: + template + DissipativeRiemannSolverExtra(FluidI& fluid_i, FluidJ& fluid_j) + : DissipativeRiemannSolver(fluid_i, fluid_j) {}; + Vecd DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij); + }; +} // namespace SPH +#endif // RIEMANN_SOLVER_EXTRA_H diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp deleted file mode 100644 index 037b7be579..0000000000 --- a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include "wetting_coupled_spatial_temporal_method.h" - -namespace SPH -{ -//=====================================================================================================// -namespace fluid_dynamics -{ -//=================================================================================================// -NonWettingSurfaceIndication:: - NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, - BaseContactRelation &contact_relation, Real threshold, Real criterion) - : FreeSurfaceIndicationComplex(inner_relation, contact_relation, threshold), wetting_criterion(criterion) -{ - for (size_t k = 0; k != contact_particles_.size(); ++k) - { - contact_phi_.push_back(this->contact_particles_[k]->template getVariableByName("Phi")); - } -} -//=================================================================================================// -NonWettingSurfaceIndication:: - NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold, Real criterion) - : NonWettingSurfaceIndication(complex_relation.getInnerRelation(), - complex_relation.getContactRelation(), threshold, criterion) {} - -//=================================================================================================// -} // namespace fluid_dynamics - //=================================================================================================// -} // namespace SPH - //=================================================================================================// \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h deleted file mode 100644 index bed6605bfd..0000000000 --- a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h +++ /dev/null @@ -1,107 +0,0 @@ -/* ------------------------------------------------------------------------- * - * 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 fluid_surface_complex.h - * @brief Here, we define the algorithm classes for fluid dynamics within the body. - * @details We consider here weakly compressible fluids. The algorithms may be - * different for free surface flow and the one without free surface. - * @author Chi Zhang and Xiangyu Hu - */ - -#ifndef WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H -#define WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H - -#include "fluid_surface_complex.h" - -namespace SPH -{ -namespace fluid_dynamics -{ -/** - * @class NonWettingSurfaceIndication - * @brief Non wetting surface particles include free-surface ones and interfacial ones near the non-wetted structure. - * @brief Even the position divergence of interfacial fluid pariticles has satisfied with the threshold of spatial-temporal - identification approach to be identified as internal ones,they will remain as free-surface ones if without - any wetted neighboring solid particles. - */ -class NonWettingSurfaceIndication : public FreeSurfaceIndicationComplex -{ - public: - NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, - BaseContactRelation &contact_relation, Real threshold = 0.75, Real criterion = 0.99); - explicit NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold = 0.75, Real criterion = 0.99); - virtual ~NonWettingSurfaceIndication(){}; - - inline void interaction(size_t index_i, Real dt = 0.0) - { - FreeSurfaceIndicationInner::interaction(index_i, dt); - - Real pos_div = 0.0; - 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) - { - pos_div -= contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.r_ij_[n]; - } - } - pos_div_[index_i] += pos_div; - - if (pos_div_[index_i] > this->threshold_by_dimensions_) - { - 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) - { - size_t j = contact_neighborhood.j_[n]; - if ((*(contact_phi_[k]))[j] > wetting_criterion) - { - pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_; - break; - } - else - { - pos_div_[index_i] = 0.5 * this->threshold_by_dimensions_; - } - } - if (pos_div_[index_i] == 2.0 * this->threshold_by_dimensions_) - break; - } - } - }; - - protected: - Real wetting_criterion; - StdVec *> contact_phi_; -}; - -using WettingCoupledSpatialTemporalFreeSurfaceIdentificationComplex = - SpatialTemporalFreeSurfaceIdentification; -using SpatialTemporalFreeSurfaceIdentificationComplex = - SpatialTemporalFreeSurfaceIdentification; - - -} // namespace fluid_dynamics -} // namespace SPH -#endif // WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H \ No newline at end of file diff --git a/tests/user_examples/test_2d_column_collapse/CMakeLists.txt b/tests/user_examples/test_2d_column_collapse/CMakeLists.txt new file mode 100644 index 0000000000..3af25fce68 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/CMakeLists.txt @@ -0,0 +1,21 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ + DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} extra_sources_2d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME} + COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/user_examples/test_2d_column_collapse/column_collapse.cpp b/tests/user_examples/test_2d_column_collapse/column_collapse.cpp new file mode 100644 index 0000000000..61196d440d --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/column_collapse.cpp @@ -0,0 +1,278 @@ +/** + * @file column_collapse.cpp + * @brief 2D dambreak example. + * @details This is the one of the basic test cases, also the first case for + * understanding SPH method for soil simulation. + */ +#include "sphinxsys.h" //SPHinXsys Library. +#include "all_continuum.h" +using namespace SPH; // Namespace cite here. +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +//unit system - 1 +Real DL = 0.5; /**< Tank length. */ +Real DH = 0.15; /**< Tank height. */ +Real LL = 0.2; /**< Liquid column length. */ +Real LH = 0.1; /**< Liquid column height. */ + +Real particle_spacing_ref = LH / 50; /**< Initial reference particle spacing. */ +Real BW = particle_spacing_ref * 4; /**< Extending width for boundary conditions. */ +BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(DL + BW, DH + BW)); +// observer location +StdVec observation_location = { Vecd(DL, 0.2) }; +//---------------------------------------------------------------------- +// Material properties of the soil. +//---------------------------------------------------------------------- +/* + * Dilatancy angle is always zero for non-associate flow rule + */ +Real rho0_s = 2650; /**< Reference density of soil. */ +Real gravity_g = 9.8; /**< Gravity force of soil. */ +Real Youngs_modulus = 0.84e6; //reference Youngs modulus +Real poisson = 0.3; //Poisson ratio +Real c_s = sqrt(Youngs_modulus / (rho0_s*3*(1-2* poisson))); +Real friction_angle = 19.8 * Pi / 180; +//Real cohesion = 100; //tensile instability occurs +Real cohesion = 0.0; +Real dilatancy = 0.0 * Pi / 180; +//---------------------------------------------------------------------- +// Geometric shapes used in this case. +//---------------------------------------------------------------------- +Vec2d soil_block_halfsize = Vec2d(0.5 * LL, 0.5 * LH); // local center at origin: +Vec2d soil_block_translation = soil_block_halfsize; +Vec2d outer_wall_halfsize = Vec2d(0.5 * DL + BW, 0.5 * DH + BW); +Vec2d outer_wall_translation = Vec2d(-BW, -BW) + outer_wall_halfsize; +Vec2d inner_wall_halfsize = Vec2d(0.5 * DL, 0.5 * DH); +Vec2d inner_wall_translation = inner_wall_halfsize; +//---------------------------------------------------------------------- +// Complex for wall boundary +//---------------------------------------------------------------------- +class WallBoundary : public ComplexShape +{ +public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + add>(Transform(outer_wall_translation), outer_wall_halfsize); + subtract>(Transform(inner_wall_translation), inner_wall_halfsize); + } +}; +std::vector soil_shape{ + Vecd(0, 0), Vecd(0, LH), Vecd(LL, LH), Vecd(LL, 0), Vecd(0, 0) }; + +class Soil : public MultiPolygonShape +{ +public: + explicit Soil(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(soil_shape, ShapeBooleanOps::add); + } +}; +//---------------------------------------------------------------------- +// application dependent initial condition +//---------------------------------------------------------------------- +class SoilInitialCondition : public continuum_dynamics::ContinuumInitialCondition +{ +public: + explicit SoilInitialCondition(RealBody& granular_column) + : continuum_dynamics::ContinuumInitialCondition(granular_column) {}; +protected: + void update(size_t index_i, Real dt) + { + /** initial stress */ + Real y = pos_[index_i][1]; + Real gama = 1 - sin(friction_angle); + Real stress_yy = - rho0_s * gravity_g * y; + stress_tensor_3D_[index_i](1, 1) = stress_yy; + stress_tensor_3D_[index_i](0, 0) = stress_yy * gama; + stress_tensor_3D_[index_i](2, 2) = stress_yy * gama; + }; +}; +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem. + //---------------------------------------------------------------------- + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); +#ifdef BOOST_AVAILABLE + // handle command line arguments + sph_system.handleCommandlineOptions(ac, av); +#endif //---------------------------------------------------------------------- + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + RealBody soil_block( + sph_system, makeShared("GranularBody")); + //soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle, cohesion, dilatancy); + soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle); + soil_block.generateParticles(); + soil_block.addBodyStateForRecording("Pressure"); + soil_block.addBodyStateForRecording("Density"); + soil_block.addBodyStateForRecording("VerticalStress"); + soil_block.addBodyStateForRecording("AccDeviatoricPlasticStrain"); + + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + wall_boundary.addBodyStateForRecording("NormalDirection"); + + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + //---------------------------------------------------------------------- + ComplexRelation soil_block_complex(soil_block, { &wall_boundary }); + //---------------------------------------------------------------------- + // Define the main numerical methods used in the simulation. + // Note that there may be data dependence on the constructors of these methods. + //---------------------------------------------------------------------- + SimpleDynamics soil_initial_condition(soil_block); + Gravity gravity(Vecd(0.0, -gravity_g)); + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + SharedPtr gravity_ptr = makeShared(Vecd(0.0, -gravity_g)); + SimpleDynamics soil_step_initialization(soil_block, gravity_ptr); + ReduceDynamics soil_acoustic_time_step(soil_block, 0.2); + InteractionWithUpdate soil_density_by_summation(soil_block_complex); + InteractionDynamics stress_diffusion(soil_block_complex.getInnerRelation()); + //stress relaxation with Riemann solver + Dynamics1Level granular_stress_relaxation_1st(soil_block_complex); + Dynamics1Level granular_stress_relaxation_2nd(soil_block_complex); + //---------------------------------------------------------------------- + // Define the methods for I/O operations, observations + // and regression tests of the simulation. + //---------------------------------------------------------------------- + IOEnvironment io_environment(sph_system); + BodyStatesRecordingToVtp body_states_recording(io_environment, sph_system.real_bodies_); + RestartIO restart_io(io_environment, sph_system.real_bodies_); + RegressionTestDynamicTimeWarping>> + write_mechanical_energy(io_environment, soil_block, gravity_ptr); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + wall_boundary_normal_direction.exec(); + //soil_initial_condition.exec(); //It's better to apply initial condition + //---------------------------------------------------------------------- + // Load restart file if necessary. + //---------------------------------------------------------------------- + if (sph_system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep()); + soil_block.updateCellLinkedList(); + soil_block_complex.updateConfiguration(); + } + //---------------------------------------------------------------------- + // Setup for time-stepping control + //---------------------------------------------------------------------- + size_t number_of_iterations = sph_system.RestartStep(); + int screen_output_interval = 500; + int observation_sample_interval = screen_output_interval * 2; + int restart_output_interval = screen_output_interval * 10; + Real End_Time = 1.0; /**< End time. */ + Real D_Time = End_Time / 50; /**< Time stamps for output of body states. */ + Real Dt = 0.1*D_Time; + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + TimeInterval interval_computing_time_step; + TimeInterval interval_computing_soil_stress_relaxation; + TimeInterval interval_updating_configuration; + TickCount time_instance; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + body_states_recording.writeToFile(); + write_mechanical_energy.writeToFile(number_of_iterations); + //---------------------------------------------------------------------- + // Main loop starts here. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < End_Time) + { + Real integration_time = 0.0; + /** Integrate time (loop) until the next output time. */ + while (integration_time < D_Time) + { + /** outer loop for dual-time criteria time-stepping. */ + time_instance = TickCount::now(); + soil_step_initialization.exec(); + + soil_density_by_summation.exec(); + interval_computing_time_step += TickCount::now() - time_instance; + + time_instance = TickCount::now(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + Real dt = soil_acoustic_time_step.exec(); + + granular_stress_relaxation_1st.exec(dt); + stress_diffusion.exec(); + granular_stress_relaxation_2nd.exec(dt); + + + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + + interval_computing_soil_stress_relaxation += TickCount::now() - time_instance; + + /** screen output, write body reduced values and restart files */ + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << std::setprecision(4) << " Time = " + << GlobalStaticVariables::physical_time_ + << std::scientific << " dt = " << dt << "\n"; + + if (number_of_iterations % observation_sample_interval == 0 && number_of_iterations != sph_system.RestartStep()) + { + write_mechanical_energy.writeToFile(number_of_iterations); + } + if (number_of_iterations % restart_output_interval == 0) + restart_io.writeToFile(number_of_iterations); + } + number_of_iterations++; + /** Update cell linked list and configuration. */ + soil_block.updateCellLinkedList(); + soil_block_complex.updateConfiguration(); + } + time_instance = TickCount::now(); + + interval_updating_configuration += TickCount::now() - time_instance; + } + TickCount t2 = TickCount::now(); + body_states_recording.writeToFile(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << std::fixed << "Total wall time for computation: " << tt.seconds() + << " seconds." << std::endl; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_time_step =" + << interval_computing_time_step.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_soil_stress_relaxation = " + << interval_computing_soil_stress_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + + //sph_system.generate_regression_data_ = true; + if (sph_system.generate_regression_data_) + { + write_mechanical_energy.generateDataBase(1.0e-3); + } + else if (sph_system.RestartStep() == 0) + { + write_mechanical_energy.testResult(); + } + + return 0; +}; \ No newline at end of file diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml new file mode 100644 index 0000000000..5eff694723 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml @@ -0,0 +1,36 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml new file mode 100644 index 0000000000..91027bf718 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml new file mode 100644 index 0000000000..3597857c59 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml new file mode 100644 index 0000000000..fc987be77c --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml new file mode 100644 index 0000000000..c600313627 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml new file mode 100644 index 0000000000..290c1f97e8 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml new file mode 100644 index 0000000000..e9f4e83311 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat new file mode 100644 index 0000000000..5cc2fa0e29 --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat @@ -0,0 +1,3 @@ +true +12 +4 \ No newline at end of file diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_column_collapse/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..56acde80de --- /dev/null +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/regression_test_tool.py @@ -0,0 +1,37 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_2d_column_collapse +""" + +case_name = "test_2d_column_collapse" +body_name = "GranularBody" +parameter_name = "TotalMechanicalEnergy" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if converged == "true": + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/user_examples/test_3d_repose_angle/CMakeLists.txt b/tests/user_examples/test_3d_repose_angle/CMakeLists.txt new file mode 100644 index 0000000000..0ee37144e2 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/CMakeLists.txt @@ -0,0 +1,21 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ + DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} extra_sources_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME} + COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/FluidObserver_Pressure_runtimes.dat b/tests/user_examples/test_3d_repose_angle/regression_test_tool/FluidObserver_Pressure_runtimes.dat new file mode 100644 index 0000000000..0db62092e1 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/FluidObserver_Pressure_runtimes.dat @@ -0,0 +1,3 @@ +false +8 +0 \ No newline at end of file diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml new file mode 100644 index 0000000000..6ef869bfd6 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml new file mode 100644 index 0000000000..443454a315 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml new file mode 100644 index 0000000000..f073ca66b0 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml new file mode 100644 index 0000000000..37dfbf4527 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml new file mode 100644 index 0000000000..74dcd4fab1 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml new file mode 100644 index 0000000000..63624ad3de --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml new file mode 100644 index 0000000000..8a87f549e6 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml new file mode 100644 index 0000000000..243776f726 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml new file mode 100644 index 0000000000..f52478f033 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat new file mode 100644 index 0000000000..34de7daedb --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat @@ -0,0 +1,3 @@ +true +7 +4 \ No newline at end of file diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_3d_repose_angle/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..99f1d04f24 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/regression_test_tool.py @@ -0,0 +1,38 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + + +""" +case name: test_3d_repose_angle +""" + +case_name = "test_3d_repose_angle" +body_name = "GranularBody" +parameter_name = "TotalMechanicalEnergy" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true"): + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/user_examples/test_3d_repose_angle/repose_angle.cpp b/tests/user_examples/test_3d_repose_angle/repose_angle.cpp new file mode 100644 index 0000000000..7735927a45 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/repose_angle.cpp @@ -0,0 +1,295 @@ +#include "sphinxsys.h" // SPHinXsys Library. +#include "all_continuum.h" +using namespace SPH; +// general parameters for geometry + +Real radius = 0.1; // liquid length +Real height = 0.1; // liquid height +Real resolution_ref = radius / 10; // particle spacing +Real BW = resolution_ref * 4; // boundary width +Real DL = 2 * radius*(1+1.24*height/radius) + 0.1; // tank length +Real DH = height + 0.02; // tank height +Real DW = DL; // tank width +// for material properties of the fluid +Real rho0_s = 2600; /**< Reference density of soil. */ +Real gravity_g = 9.8; /**< Gravity force of soil. */ +Real Youngs_modulus = 5.98e6; //reference Youngs modulus +Real poisson = 0.3; //Poisson ratio +Real c_s = sqrt(Youngs_modulus / (rho0_s * 3 * (1 - 2 * poisson))); +Real friction_angle = 30 * Pi / 180; +Real cohesion = 0; +Real dilatancy = 0 * Pi / 180; +/** Define the soil body. */ +Real inner_circle_radius = radius; +int resolution(20); +class SoilBlock : public ComplexShape +{ +public: + explicit SoilBlock(const std::string& shape_name) : ComplexShape(shape_name) + { + Vecd translation_column(DL / 2, 0.5 * height, DW / 2); + add(SimTK::UnitVec3(0, 1.0, 0), inner_circle_radius, + 0.5 * height, resolution, translation_column); + } +}; +// define the static solid wall boundary shape +class WallBoundary : public ComplexShape +{ +public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + Vecd outer_wall_halfsize = Vecd(0.5 * DL + BW, 0.5 * DH + BW, 0.5 * DW + BW); + Vecd outer_wall_translation = Vecd(-BW, -BW, -BW) + outer_wall_halfsize; + Vecd inner_wall_halfsize = Vecd(0.5 * DL, 0.5 * DH, 0.5 * DW); + Vecd inner_wall_translation = inner_wall_halfsize; + add>(Transform(outer_wall_translation), outer_wall_halfsize); + subtract>(Transform(inner_wall_translation), inner_wall_halfsize); + } +}; +//---------------------------------------------------------------------- +// application dependent initial condition +//---------------------------------------------------------------------- +class SoilInitialCondition : public continuum_dynamics::ContinuumInitialCondition +{ +public: + explicit SoilInitialCondition(RealBody& granular_column) + : continuum_dynamics::ContinuumInitialCondition(granular_column) {}; +protected: + void update(size_t index_i, Real dt) + { + /** initial stress */ + Real y = pos_[index_i][1]; + //Real gama = poisson / (1 - poisson); + Real gama = 1 - sin(friction_angle); + Real stress_yy = -rho0_s * gravity_g * y; + stress_tensor_3D_[index_i](1, 1) = stress_yy; + stress_tensor_3D_[index_i](0, 0) = stress_yy * gama; + stress_tensor_3D_[index_i](2, 2) = stress_yy * gama; + }; +}; +// the main program with commandline options +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up an SPHSystem. + //---------------------------------------------------------------------- + BoundingBox system_domain_bounds(Vecd(-BW, -BW, -BW), Vecd(DL + BW, DH + BW, DW + BW)); + SPHSystem sph_system(system_domain_bounds, resolution_ref); + sph_system.setRunParticleRelaxation(false); + sph_system.setReloadParticles(true); + sph_system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(sph_system); + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + RealBody soil_block(sph_system, makeShared("GranularBody")); + soil_block.defineBodyLevelSetShape()->writeLevelSet(io_environment); + soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle, cohesion, dilatancy); + (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) + ? soil_block.generateParticles(io_environment, soil_block.getName()) + : soil_block.generateParticles(); + soil_block.addBodyStateForRecording("Pressure"); + soil_block.addBodyStateForRecording("Density"); + soil_block.addBodyStateForRecording("VerticalStress"); + soil_block.addBodyStateForRecording("AccDeviatoricPlasticStrain"); + + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + wall_boundary.addBodyStateForRecording("NormalDirection"); + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + //---------------------------------------------------------------------- + ComplexRelation soil_block_complex(soil_block, { &wall_boundary }); + BodyStatesRecordingToVtp body_states_recording(io_environment, sph_system.real_bodies_); + //run particle relaxation + if (sph_system.RunParticleRelaxation()) + { + /** + * @brief Methods used for particle relaxation. + */ + /** Random reset the insert body particle position. */ + SimpleDynamics random_column_particles(soil_block); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_column_to_vtp(io_environment, soil_block); + /** Write the particle reload files. */ + + ReloadParticleIO write_particle_reload_files(io_environment, soil_block); + /** A Physics relaxation step. */ + relax_dynamics::RelaxationStepInner relaxation_step_inner(soil_block_complex.getInnerRelation()); + /** + * @brief Particle relaxation starts here. + */ + random_column_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + body_states_recording.writeToFile(0.0); + + /** relax particles of the insert body. */ + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; + write_column_to_vtp.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of cylinder body finish !" << std::endl; + /** Output results. */ + write_particle_reload_files.writeToFile(0.0); + return 0; + } + //---------------------------------------------------------------------- + // Define the numerical methods used in the simulation. + // Note that there may be data dependence on the sequence of constructions. + //---------------------------------------------------------------------- + SimpleDynamics soil_initial_condition(soil_block); + SharedPtr gravity_ptr = makeShared(Vec3d(0.0, -gravity_g, 0.0)); + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + SimpleDynamics soil_step_initialization(soil_block, gravity_ptr); + ReduceDynamics soil_acoustic_time_step(soil_block, 0.1); + InteractionWithUpdate soil_density_by_summation(soil_block_complex); + InteractionDynamics stress_diffusion(soil_block_complex.getInnerRelation()); + Dynamics1Level granular_stress_relaxation_1st(soil_block_complex); + Dynamics1Level granular_stress_relaxation_2nd(soil_block_complex); + //---------------------------------------------------------------------- + // Define the methods for I/O operations, observations + // and regression tests of the simulation. + //---------------------------------------------------------------------- + RestartIO restart_io(io_environment, sph_system.real_bodies_); + RegressionTestDynamicTimeWarping>> + write_soil_mechanical_energy(io_environment, soil_block, gravity_ptr); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + wall_boundary_normal_direction.exec(); + soil_initial_condition.exec(); + //---------------------------------------------------------------------- + // Load restart file if necessary. + //---------------------------------------------------------------------- + if (sph_system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep()); + soil_block.updateCellLinkedList(); + soil_block_complex.updateConfiguration(); + } + //---------------------------------------------------------------------- + // Setup for time-stepping control + //---------------------------------------------------------------------- + size_t number_of_iterations = sph_system.RestartStep(); + int screen_output_interval = 500; + int observation_sample_interval = screen_output_interval * 2; + int restart_output_interval = screen_output_interval * 10; + Real End_Time = 0.5; /**< End time. */ + Real D_Time = End_Time / 25; /**< Time stamps for output of body states. */ + Real Dt = 0.1 * D_Time; + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + TimeInterval interval_computing_time_step; + TimeInterval interval_computing_soil_stress_relaxation; + TimeInterval interval_updating_configuration; + TickCount time_instance; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + body_states_recording.writeToFile(); + write_soil_mechanical_energy.writeToFile(number_of_iterations); + //---------------------------------------------------------------------- + // Main loop starts here. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < End_Time) + { + Real integration_time = 0.0; + /** Integrate time (loop) until the next output time. */ + while (integration_time < D_Time) + { + /** outer loop for dual-time criteria time-stepping. */ + time_instance = TickCount::now(); + soil_step_initialization.exec(); + + soil_density_by_summation.exec(); + interval_computing_time_step += TickCount::now() - time_instance; + + time_instance = TickCount::now(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + Real dt = soil_acoustic_time_step.exec(); + + granular_stress_relaxation_1st.exec(dt); + stress_diffusion.exec(); + granular_stress_relaxation_2nd.exec(dt); + + + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + + interval_computing_soil_stress_relaxation += TickCount::now() - time_instance; + + /** screen output, write body reduced values and restart files */ + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << std::setprecision(4) << " Time = " + << GlobalStaticVariables::physical_time_ + << std::scientific << " dt = " << dt << "\n"; + + if (number_of_iterations % observation_sample_interval == 0 && number_of_iterations != sph_system.RestartStep()) + { + write_soil_mechanical_energy.writeToFile(number_of_iterations); + } + if (number_of_iterations % restart_output_interval == 0) + restart_io.writeToFile(number_of_iterations); + } + number_of_iterations++; + + soil_block.updateCellLinkedList(); + soil_block_complex.updateConfiguration(); + } + /** Update cell linked list and configuration. */ + time_instance = TickCount::now(); + interval_updating_configuration += TickCount::now() - time_instance; + } + + TickCount t2 = TickCount::now(); + body_states_recording.writeToFile(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << std::fixed << "Total wall time for computation: " << tt.seconds() + << " seconds." << std::endl; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_time_step =" + << interval_computing_time_step.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_soil_stress_relaxation = " + << interval_computing_soil_stress_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + std::cout << "total time steps = " << number_of_iterations << "\n"; + + //sph_system.generate_regression_data_ = true; + if (sph_system.generate_regression_data_) + { + write_soil_mechanical_energy.generateDataBase(1.0e-3); + } + else if (sph_system.RestartStep() == 0) + { + write_soil_mechanical_energy.testResult(); + } + + return 0; +} From 7d3e0b40415ff41b261ccc9fcf274a4e4d9bbea9 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 11:52:37 +0200 Subject: [PATCH 2/7] clean the code and add regression test --- .../continuum_dynamics_complex.hpp | 15 ----------- .../continuum_dynamics_inner.cpp | 7 ----- .../continuum_dynamics_inner.h | 4 --- .../continuum_dynamics_inner.hpp | 16 ----------- .../general_continuum/continuum_particles.cpp | 3 --- .../column_collapse.cpp | 27 ------------------- ...Body_TotalMechanicalEnergy_dtwdistance.xml | 2 +- .../test_3d_repose_angle/repose_angle.cpp | 13 ++------- 8 files changed, 3 insertions(+), 84 deletions(-) diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp index 6f86f35f56..43ac0596cc 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp @@ -87,14 +87,10 @@ void BaseStressRelaxation1stHalfWithWall::inter BaseStressRelaxation1stHalfType::interaction(index_i, dt); Vecd acc_prior_i = computeNonConservativeAcceleration(index_i); - - //Vecd acceleration = Vecd::Zero(); Vecd acceleration = acc_prior_i; Real rho_dissipation(0); Matd stress_tensor_i = this->reduceTensor(this->stress_tensor_3D_[index_i]); - //Real rho_i = this->rho_[index_i]; - //Real rho_in_wall = this->plastic_continuum_.getDensity(); for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) { @@ -110,15 +106,11 @@ void BaseStressRelaxation1stHalfWithWall::inter Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_wall_external_acceleration); - //Real p_in_wall = this->p_[index_i]; Matd stress_tensor_in_wall = stress_tensor_i; - //acceleration += rho_in_wall * (stress_tensor_i / (rho_i * rho_i) + stress_tensor_in_wall / (rho_in_wall * rho_in_wall)) * nablaW_ijV_j; acceleration += (stress_tensor_i + stress_tensor_in_wall) * nablaW_ijV_j; - rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; } } - //this->acc_[index_i] += acceleration; this->acc_[index_i] += acceleration / this->rho_[index_i]; this->drho_dt_[index_i] += rho_dissipation * this->rho_[index_i]; } @@ -154,19 +146,14 @@ void BaseStressRelaxation2ndHalfWithWall::inter velocity_gradient += velocity_gradient_ij; density_change_rate += (this->vel_[index_i] - vel_in_wall).dot(e_ij) * dW_ijV_j; - // Real u_jump = 2.0 * (this->vel_[index_i] - vel_ave_k[index_j]).dot(n_k[index_j]); - // p_dissipation += this->riemann_solver_.DissipativePJump(u_jump) * dW_ijV_j * n_k[index_j]; Vecd u_jump = 2.0 * (this->vel_[index_i] - vel_ave_k[index_j]); p_dissipation += this->riemann_solver_.DissipativePJumpExtra(u_jump, n_k[index_j]) * dW_ijV_j; - // Real temp = u_jump.dot(n_k[index_j]); - // p_dissipation += this->riemann_solver_.DissipativePJump(temp) * dW_ijV_j * e_ij; } } this->drho_dt_[index_i] += density_change_rate * this->rho_[index_i]; this->velocity_gradient_[index_i] += velocity_gradient; this->acc_[index_i] += p_dissipation / this->rho_[index_i]; } - //=================================================================================================// //================================BaseStressDiffusionWithWall======================================// //=================================================================================================// @@ -193,12 +180,10 @@ void BaseStressDiffusionWithWall::interaction(size_t in Real y_ij = this->pos_[index_i](1, 0) - this->pos_[index_j](1, 0); //stress boundary condition Mat3d stress_tensor_j = stress_tensor_i; - diffusion_stress_ = stress_tensor_i - stress_tensor_j; diffusion_stress_(0, 0) = diffusion_stress_(0, 0) - (1 - sin(this->fai_)) * density * gravity * y_ij; diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij; diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(this->fai_)) * density * gravity * y_ij; - diffusion_stress_rate_ += 2 * this->zeta_ * this->smoothing_length_ * this->sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * this->smoothing_length_); } } diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp index 39a96314ed..6302251f51 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp @@ -245,7 +245,6 @@ namespace SPH Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); } - //=================================================================================================// FixedInAxisDirection::FixedInAxisDirection(BodyPartByParticle& body_part, Vecd constrained_axises) : BaseMotionConstraint(body_part), constrain_matrix_(Matd::Identity()) @@ -258,7 +257,6 @@ namespace SPH { vel_[index_i] = constrain_matrix_ * vel_[index_i]; }; - //=================================================================================================// ConstrainSolidBodyMassCenter:: ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction) @@ -304,7 +302,6 @@ namespace SPH } return tensor_2d; } - Mat3d BaseRelaxationPlastic::increaseTensor(Matd tensor_2d) { Mat3d tensor_3d = Mat3d::Zero(); @@ -317,7 +314,6 @@ namespace SPH } return tensor_3d; } - //====================================================================================// //===============================StressDiffusion======================================// //====================================================================================// @@ -338,13 +334,10 @@ namespace SPH Real r_ij = inner_neighborhood.r_ij_[n]; Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); - //calculate Psi, equation (34) in Feng_2021_CG diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; diffusion_stress_(0, 0) = diffusion_stress_(0, 0) - (1 - sin(fai_)) * density * gravity * y_ij; diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij; diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(fai_)) * density * gravity * y_ij; - - //diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_ * smoothing_length_); diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); } stress_rate_3D_[index_i] = diffusion_stress_rate_; diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h index 2ccc7d8714..dbad3d1932 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h @@ -164,8 +164,6 @@ class AngularConservativeShearAccelerationRelaxation : public ShearAccelerationR virtual ~AngularConservativeShearAccelerationRelaxation() {}; void interaction(size_t index_i, Real dt = 0.0); - //void update(size_t index_i, Real dt = 0.0); - }; /** @@ -278,7 +276,6 @@ class BaseRelaxationPlastic : public LocalDynamics, public PlasticContinuumDataI StdLargeVec& pos_, & vel_, & acc_, & acc_prior_; StdLargeVec& stress_tensor_3D_, & strain_tensor_3D_, & stress_rate_3D_, & strain_rate_3D_; StdLargeVec& elastic_strain_tensor_3D_, & elastic_strain_rate_3D_; - }; //=================================================================================================// @@ -333,7 +330,6 @@ class BaseStressRelaxation2ndHalf : public BaseRelaxationPlastic StdLargeVec& Vol_, & mass_, & von_mises_stress_; StdLargeVec& acc_deviatoric_plastic_strain_, & vertical_stress_; Real E_, nu_; - }; using StressRelaxation2ndHalf = BaseStressRelaxation2ndHalf; using StressRelaxation2ndHalfRiemann = BaseStressRelaxation2ndHalf; diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp index cb626b7f93..028d980a00 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp @@ -38,15 +38,9 @@ namespace SPH { rho_[index_i] += drho_dt_[index_i] * dt * 0.5; p_[index_i] = plastic_continuum_.getPressure(rho_[index_i]); - //pressure_replace_by_volumn_part - //p_[index_i] = - stress_tensor_3D_[index_i].trace() /3; - pos_[index_i] += vel_[index_i] * dt * 0.5; - stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt * 0.5; strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt * 0.5; - - //calculate elastic strain elastic_strain_tensor_3D_[index_i] += elastic_strain_rate_3D_[index_i] * dt * 0.5; } //=================================================================================================// @@ -87,7 +81,6 @@ namespace SPH Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; Matd stress_tensor_j = reduceTensor(stress_tensor_3D_[index_j]); - //acceleration += rho_[index_j] * (stress_tensor_i / (rho_i * rho_i) + stress_tensor_j / (rho_[index_j] * rho_[index_j])) * nablaW_ijV_j; acceleration += rho_[index_j] * ((stress_tensor_i + stress_tensor_j) / (rho_i * rho_[index_j])) * nablaW_ijV_j; rho_dissipation += riemann_solver_.DissipativeUJump(p_[index_i] - p_[index_j]) * dW_ijV_j; } @@ -123,9 +116,6 @@ namespace SPH Vecd u_jump = vel_[index_i] - vel_[index_j]; density_change_rate += u_jump.dot(e_ij) * dW_ijV_j; p_dissipation += riemann_solver_.DissipativePJumpExtra(u_jump, e_ij) * dW_ijV_j; - // Real temp = u_jump.dot(e_ij); - // p_dissipation += riemann_solver_.DissipativePJump(temp) * dW_ijV_j * e_ij; - Matd velocity_gradient_ij = -(vel_[index_i] - vel_[index_j]) * nablaW_ijV_j.transpose(); velocity_gradient += velocity_gradient_ij; } @@ -142,18 +132,12 @@ namespace SPH Mat3d velocity_gradient = increaseTensor(velocity_gradient_[index_i]); Mat3d stress_tensor_rate_3D_ = plastic_continuum_.ConstitutiveRelationZ(velocity_gradient, stress_tensor_3D_[index_i]); - //consider diffusion stress with += - //stress_rate_3D_[index_i] = stress_tensor_rate_3D_; stress_rate_3D_[index_i] += stress_tensor_rate_3D_; - //update stress tensor stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt * 0.5; - //For plasticity Mat3d stress_tensor_ = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); stress_tensor_3D_[index_i] = stress_tensor_; - //added vertical_stress_[index_i] = stress_tensor_3D_[index_i](1, 1); - //Accumulated deviatoric plastic strains strain_rate_3D_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt * 0.5; //calculate elastic strain diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp index 4c169dd08c..463ebfde57 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp @@ -23,7 +23,6 @@ namespace SPH //---------------------------------------------------------------------- // register sortable particle data //---------------------------------------------------------------------- - registerSortableVariable("AccelerationByShear"); registerSortableVariable("StressTensor"); registerSortableVariable("StressTensorRate"); @@ -42,7 +41,6 @@ namespace SPH [&](size_t i) -> Vecd { return n_[i]; }); } - //=================================================================================================// PlasticContinuumParticles:: PlasticContinuumParticles(SPHBody& sph_body, PlasticContinuum* plastic_continuum) @@ -81,7 +79,6 @@ namespace SPH registerSortableVariable("VerticalStress"); registerSortableVariable("AccDeviatoricPlasticStrain"); - } //=================================================================================================// Real PlasticContinuumParticles::getDeviatoricPlasticStrain(Mat3d& strain_tensor) diff --git a/tests/user_examples/test_2d_column_collapse/column_collapse.cpp b/tests/user_examples/test_2d_column_collapse/column_collapse.cpp index 61196d440d..b6284133fc 100644 --- a/tests/user_examples/test_2d_column_collapse/column_collapse.cpp +++ b/tests/user_examples/test_2d_column_collapse/column_collapse.cpp @@ -33,9 +33,6 @@ Real Youngs_modulus = 0.84e6; //reference Youngs modulus Real poisson = 0.3; //Poisson ratio Real c_s = sqrt(Youngs_modulus / (rho0_s*3*(1-2* poisson))); Real friction_angle = 19.8 * Pi / 180; -//Real cohesion = 100; //tensile instability occurs -Real cohesion = 0.0; -Real dilatancy = 0.0 * Pi / 180; //---------------------------------------------------------------------- // Geometric shapes used in this case. //---------------------------------------------------------------------- @@ -69,26 +66,6 @@ class Soil : public MultiPolygonShape } }; //---------------------------------------------------------------------- -// application dependent initial condition -//---------------------------------------------------------------------- -class SoilInitialCondition : public continuum_dynamics::ContinuumInitialCondition -{ -public: - explicit SoilInitialCondition(RealBody& granular_column) - : continuum_dynamics::ContinuumInitialCondition(granular_column) {}; -protected: - void update(size_t index_i, Real dt) - { - /** initial stress */ - Real y = pos_[index_i][1]; - Real gama = 1 - sin(friction_angle); - Real stress_yy = - rho0_s * gravity_g * y; - stress_tensor_3D_[index_i](1, 1) = stress_yy; - stress_tensor_3D_[index_i](0, 0) = stress_yy * gama; - stress_tensor_3D_[index_i](2, 2) = stress_yy * gama; - }; -}; -//---------------------------------------------------------------------- // Main program starts here. //---------------------------------------------------------------------- int main(int ac, char *av[]) @@ -106,7 +83,6 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- RealBody soil_block( sph_system, makeShared("GranularBody")); - //soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle, cohesion, dilatancy); soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle); soil_block.generateParticles(); soil_block.addBodyStateForRecording("Pressure"); @@ -129,7 +105,6 @@ int main(int ac, char *av[]) // Define the main numerical methods used in the simulation. // Note that there may be data dependence on the constructors of these methods. //---------------------------------------------------------------------- - SimpleDynamics soil_initial_condition(soil_block); Gravity gravity(Vecd(0.0, -gravity_g)); SimpleDynamics wall_boundary_normal_direction(wall_boundary); SharedPtr gravity_ptr = makeShared(Vecd(0.0, -gravity_g)); @@ -156,7 +131,6 @@ int main(int ac, char *av[]) sph_system.initializeSystemCellLinkedLists(); sph_system.initializeSystemConfigurations(); wall_boundary_normal_direction.exec(); - //soil_initial_condition.exec(); //It's better to apply initial condition //---------------------------------------------------------------------- // Load restart file if necessary. //---------------------------------------------------------------------- @@ -216,7 +190,6 @@ int main(int ac, char *av[]) stress_diffusion.exec(); granular_stress_relaxation_2nd.exec(dt); - relaxation_time += dt; integration_time += dt; GlobalStaticVariables::physical_time_ += dt; diff --git a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml index f52478f033..04adca364f 100644 --- a/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml +++ b/tests/user_examples/test_3d_repose_angle/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/user_examples/test_3d_repose_angle/repose_angle.cpp b/tests/user_examples/test_3d_repose_angle/repose_angle.cpp index 7735927a45..21901e4d99 100644 --- a/tests/user_examples/test_3d_repose_angle/repose_angle.cpp +++ b/tests/user_examples/test_3d_repose_angle/repose_angle.cpp @@ -2,7 +2,6 @@ #include "all_continuum.h" using namespace SPH; // general parameters for geometry - Real radius = 0.1; // liquid length Real height = 0.1; // liquid height Real resolution_ref = radius / 10; // particle spacing @@ -10,15 +9,13 @@ Real BW = resolution_ref * 4; // boundary width Real DL = 2 * radius*(1+1.24*height/radius) + 0.1; // tank length Real DH = height + 0.02; // tank height Real DW = DL; // tank width -// for material properties of the fluid +// for material properties Real rho0_s = 2600; /**< Reference density of soil. */ Real gravity_g = 9.8; /**< Gravity force of soil. */ Real Youngs_modulus = 5.98e6; //reference Youngs modulus Real poisson = 0.3; //Poisson ratio Real c_s = sqrt(Youngs_modulus / (rho0_s * 3 * (1 - 2 * poisson))); Real friction_angle = 30 * Pi / 180; -Real cohesion = 0; -Real dilatancy = 0 * Pi / 180; /** Define the soil body. */ Real inner_circle_radius = radius; int resolution(20); @@ -59,7 +56,6 @@ class SoilInitialCondition : public continuum_dynamics::ContinuumInitialConditio { /** initial stress */ Real y = pos_[index_i][1]; - //Real gama = poisson / (1 - poisson); Real gama = 1 - sin(friction_angle); Real stress_yy = -rho0_s * gravity_g * y; stress_tensor_3D_[index_i](1, 1) = stress_yy; @@ -84,7 +80,7 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- RealBody soil_block(sph_system, makeShared("GranularBody")); soil_block.defineBodyLevelSetShape()->writeLevelSet(io_environment); - soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle, cohesion, dilatancy); + soil_block.defineParticlesAndMaterial(rho0_s, c_s, Youngs_modulus, poisson, friction_angle); (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) ? soil_block.generateParticles(io_environment, soil_block.getName()) : soil_block.generateParticles(); @@ -230,7 +226,6 @@ int main(int ac, char *av[]) stress_diffusion.exec(); granular_stress_relaxation_2nd.exec(dt); - relaxation_time += dt; integration_time += dt; GlobalStaticVariables::physical_time_ += dt; @@ -252,7 +247,6 @@ int main(int ac, char *av[]) restart_io.writeToFile(number_of_iterations); } number_of_iterations++; - soil_block.updateCellLinkedList(); soil_block_complex.updateConfiguration(); } @@ -260,12 +254,10 @@ int main(int ac, char *av[]) time_instance = TickCount::now(); interval_updating_configuration += TickCount::now() - time_instance; } - TickCount t2 = TickCount::now(); body_states_recording.writeToFile(); TickCount t3 = TickCount::now(); interval += t3 - t2; - } TickCount t4 = TickCount::now(); @@ -290,6 +282,5 @@ int main(int ac, char *av[]) { write_soil_mechanical_energy.testResult(); } - return 0; } From cae7407a0947b2eb3a3079c0577bc9ae5213d3b8 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 12:24:07 +0200 Subject: [PATCH 3/7] fix a bug on MacOS --- .../continuum_dynamics/continuum_dynamics_inner.hpp | 2 +- .../shared/general_continuum/general_continuum.cpp | 6 +----- .../extra_src/shared/general_continuum/general_continuum.h | 3 +-- 3 files changed, 3 insertions(+), 8 deletions(-) diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp index 028d980a00..43556b184e 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp @@ -131,7 +131,7 @@ namespace SPH Vol_[index_i] = mass_[index_i] / rho_[index_i]; Mat3d velocity_gradient = increaseTensor(velocity_gradient_[index_i]); - Mat3d stress_tensor_rate_3D_ = plastic_continuum_.ConstitutiveRelationZ(velocity_gradient, stress_tensor_3D_[index_i]); + Mat3d stress_tensor_rate_3D_ = plastic_continuum_.ConstitutiveRelation(velocity_gradient, stress_tensor_3D_[index_i]); stress_rate_3D_[index_i] += stress_tensor_rate_3D_; stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt * 0.5; //For plasticity diff --git a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp index d247123aaa..c9f4dc0ef5 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp @@ -29,10 +29,6 @@ namespace SPH //=============================================================================================// //==============================PlasticContinuum===============================================// //=============================================================================================// - Real PlasticContinuum::getPressure(Real rho) - { - return p0_ * (rho / rho0_ - 1.0); - } Real PlasticContinuum::getDPConstantsA(Real friction_angle) { @@ -44,7 +40,7 @@ namespace SPH return 3 * cohesion / sqrt(9 + 12 * tan(friction_angle) * tan(friction_angle)); } - Mat3d PlasticContinuum::ConstitutiveRelationZ(Mat3d& velocity_gradient, Mat3d& stress_tensor) + Mat3d PlasticContinuum::ConstitutiveRelation(Mat3d& velocity_gradient, Mat3d& stress_tensor) { Real dim = 3; Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); diff --git a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h index 853c4d2442..fb31a72454 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h +++ b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h @@ -65,9 +65,8 @@ class PlasticContinuum : public GeneralContinuum Real getDPConstantsA(Real friction_angle); Real getDPConstantsK(Real cohesion, Real friction_angle); Real getFrictionAngle() { return fai_; }; - virtual Real getPressure(Real rho); - virtual Mat3d ConstitutiveRelationZ(Mat3d& velocity_gradient, Mat3d& stress_tensor); + virtual Mat3d ConstitutiveRelation(Mat3d& velocity_gradient, Mat3d& stress_tensor); virtual Mat3d ReturnMapping(Mat3d& stress_tensor); virtual GeneralContinuum* ThisObjectPtr() override { return this; }; From 3777f10fb937db5f7522b5b6b16e0d8911cf3f41 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 12:56:30 +0200 Subject: [PATCH 4/7] accidently deleted a source file from another case, already fixed it. --- ...etting_coupled_spatial_temporal_method.cpp | 29 +++++ .../wetting_coupled_spatial_temporal_method.h | 107 ++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp create mode 100644 tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp new file mode 100644 index 0000000000..037b7be579 --- /dev/null +++ b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp @@ -0,0 +1,29 @@ +#include "wetting_coupled_spatial_temporal_method.h" + +namespace SPH +{ +//=====================================================================================================// +namespace fluid_dynamics +{ +//=================================================================================================// +NonWettingSurfaceIndication:: + NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, + BaseContactRelation &contact_relation, Real threshold, Real criterion) + : FreeSurfaceIndicationComplex(inner_relation, contact_relation, threshold), wetting_criterion(criterion) +{ + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_phi_.push_back(this->contact_particles_[k]->template getVariableByName("Phi")); + } +} +//=================================================================================================// +NonWettingSurfaceIndication:: + NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold, Real criterion) + : NonWettingSurfaceIndication(complex_relation.getInnerRelation(), + complex_relation.getContactRelation(), threshold, criterion) {} + +//=================================================================================================// +} // namespace fluid_dynamics + //=================================================================================================// +} // namespace SPH + //=================================================================================================// \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h new file mode 100644 index 0000000000..bed6605bfd --- /dev/null +++ b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h @@ -0,0 +1,107 @@ +/* ------------------------------------------------------------------------- * + * 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 fluid_surface_complex.h + * @brief Here, we define the algorithm classes for fluid dynamics within the body. + * @details We consider here weakly compressible fluids. The algorithms may be + * different for free surface flow and the one without free surface. + * @author Chi Zhang and Xiangyu Hu + */ + +#ifndef WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H +#define WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H + +#include "fluid_surface_complex.h" + +namespace SPH +{ +namespace fluid_dynamics +{ +/** + * @class NonWettingSurfaceIndication + * @brief Non wetting surface particles include free-surface ones and interfacial ones near the non-wetted structure. + * @brief Even the position divergence of interfacial fluid pariticles has satisfied with the threshold of spatial-temporal + identification approach to be identified as internal ones,they will remain as free-surface ones if without + any wetted neighboring solid particles. + */ +class NonWettingSurfaceIndication : public FreeSurfaceIndicationComplex +{ + public: + NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, + BaseContactRelation &contact_relation, Real threshold = 0.75, Real criterion = 0.99); + explicit NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold = 0.75, Real criterion = 0.99); + virtual ~NonWettingSurfaceIndication(){}; + + inline void interaction(size_t index_i, Real dt = 0.0) + { + FreeSurfaceIndicationInner::interaction(index_i, dt); + + Real pos_div = 0.0; + 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) + { + pos_div -= contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.r_ij_[n]; + } + } + pos_div_[index_i] += pos_div; + + if (pos_div_[index_i] > this->threshold_by_dimensions_) + { + 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) + { + size_t j = contact_neighborhood.j_[n]; + if ((*(contact_phi_[k]))[j] > wetting_criterion) + { + pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_; + break; + } + else + { + pos_div_[index_i] = 0.5 * this->threshold_by_dimensions_; + } + } + if (pos_div_[index_i] == 2.0 * this->threshold_by_dimensions_) + break; + } + } + }; + + protected: + Real wetting_criterion; + StdVec *> contact_phi_; +}; + +using WettingCoupledSpatialTemporalFreeSurfaceIdentificationComplex = + SpatialTemporalFreeSurfaceIdentification; +using SpatialTemporalFreeSurfaceIdentificationComplex = + SpatialTemporalFreeSurfaceIdentification; + + +} // namespace fluid_dynamics +} // namespace SPH +#endif // WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H \ No newline at end of file From cc5c89f68e33e9530b6ac56402656c3fef452969 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 16:01:50 +0200 Subject: [PATCH 5/7] fix a bug when running on Linux-float --- .../continuum_dynamics_complex.hpp | 2 +- .../general_continuum/riemann_solver_extra.cpp | 3 +-- .../test_3d_repose_angle/CMakeLists.txt | 14 +++++++++++--- .../user_examples/test_3d_repose_angle/run_test.sh | 2 ++ 4 files changed, 15 insertions(+), 6 deletions(-) create mode 100644 tests/user_examples/test_3d_repose_angle/run_test.sh diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp index 43ac0596cc..1749f6d746 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp @@ -105,7 +105,7 @@ void BaseStressRelaxation1stHalfWithWall::inter Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); - Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_wall_external_acceleration); + Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(Real(0), face_wall_external_acceleration); Matd stress_tensor_in_wall = stress_tensor_i; acceleration += (stress_tensor_i + stress_tensor_in_wall) * nablaW_ijV_j; rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; diff --git a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp index 06f6f6cba7..5a6c6c87db 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp @@ -4,9 +4,8 @@ namespace SPH { Vecd AcousticRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) { - return rho0c0_geo_ave_ * u_jump * SMIN(3.0 * SMAX(u_jump.dot(e_ij) * inv_c_ave_, 0.0), 1.0); + return rho0c0_geo_ave_ * u_jump * SMIN(3.0 * SMAX(u_jump.dot(e_ij) * inv_c_ave_, Real(0)), Real(1)); } - Vecd DissipativeRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) { return rho0c0_geo_ave_ * u_jump; diff --git a/tests/user_examples/test_3d_repose_angle/CMakeLists.txt b/tests/user_examples/test_3d_repose_angle/CMakeLists.txt index 0ee37144e2..fb50cb12cf 100644 --- a/tests/user_examples/test_3d_repose_angle/CMakeLists.txt +++ b/tests/user_examples/test_3d_repose_angle/CMakeLists.txt @@ -16,6 +16,14 @@ target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) target_link_libraries(${PROJECT_NAME} extra_sources_3d) set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") -add_test(NAME ${PROJECT_NAME} - COMMAND ${PROJECT_NAME} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file +if(${CMAKE_SYSTEM_NAME} MATCHES "Windows") + add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --r=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --r=false --i=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +else() + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/run_test.sh + DESTINATION ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND bash ${EXECUTABLE_OUTPUT_PATH}/run_test.sh + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +endif() \ No newline at end of file diff --git a/tests/user_examples/test_3d_repose_angle/run_test.sh b/tests/user_examples/test_3d_repose_angle/run_test.sh new file mode 100644 index 0000000000..945e1e1ac2 --- /dev/null +++ b/tests/user_examples/test_3d_repose_angle/run_test.sh @@ -0,0 +1,2 @@ +./test_3d_repose_angle --r=true +./test_3d_repose_angle --r=false --i=true \ No newline at end of file From cfdc472f24af828cc15ba96b4ecbebc02eee65ff Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Sat, 5 Aug 2023 18:04:08 +0200 Subject: [PATCH 6/7] fix another bug on Linux-float --- .../extra_src/shared/general_continuum/riemann_solver_extra.cpp | 2 +- .../GranularBody_TotalMechanicalEnergy_dtwdistance.xml | 2 +- .../PlateBody_TotalMechanicalEnergy_dtwdistance.xml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp index 5a6c6c87db..fd8bff9aed 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/riemann_solver_extra.cpp @@ -4,7 +4,7 @@ namespace SPH { Vecd AcousticRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) { - return rho0c0_geo_ave_ * u_jump * SMIN(3.0 * SMAX(u_jump.dot(e_ij) * inv_c_ave_, Real(0)), Real(1)); + return rho0c0_geo_ave_ * u_jump * SMIN(Real(3) * SMAX(u_jump.dot(e_ij) * inv_c_ave_, Real(0)), Real(1)); } Vecd DissipativeRiemannSolverExtra::DissipativePJumpExtra(const Vecd& u_jump, const Vecd& e_ij) { diff --git a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml index e9f4e83311..22fa5165a2 100644 --- a/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml +++ b/tests/user_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalMechanicalEnergy_dtwdistance.xml index dc52c02fc2..ec66d0694c 100644 --- a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalMechanicalEnergy_dtwdistance.xml +++ b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalMechanicalEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + From b461fee770104db8cae568493d31d06e1ac6f144 Mon Sep 17 00:00:00 2001 From: RussellSyne <1477856817@qq.com> Date: Mon, 7 Aug 2023 11:40:14 +0200 Subject: [PATCH 7/7] remove redundant variables and comments --- .../continuum_dynamics_complex.h | 7 ------- .../continuum_dynamics_complex.hpp | 17 +++-------------- 2 files changed, 3 insertions(+), 21 deletions(-) diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h index c0e82470ab..e9e14ffb2a 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h @@ -57,10 +57,6 @@ class BaseShearStressRelaxation2ndHalfWithWall : public fluid_dynamics::Interact using ShearStressRelaxation2ndHalfWithWall = BaseShearStressRelaxation2ndHalfWithWall; -//=============================================Plasticity==========================================// -//=================================================================================================// -//==============================BaseStressRelaxation1stHalfWithWall=================================// -//=================================================================================================// template class BaseStressRelaxation1stHalfWithWall : public fluid_dynamics::InteractionWithWall { @@ -77,9 +73,6 @@ class BaseStressRelaxation1stHalfWithWall : public fluid_dynamics::InteractionWi using StressRelaxation1stHalfWithWall = BaseStressRelaxation1stHalfWithWall; using StressRelaxation1stHalfRiemannWithWall = BaseStressRelaxation1stHalfWithWall; using StressRelaxation1stHalfDissipativeRiemannfWithWall = BaseStressRelaxation1stHalfWithWall; -//=================================================================================================// -//===============================BaseStressRelaxation2ndHalfWithWall======================================// -//=================================================================================================// template class BaseStressRelaxation2ndHalfWithWall : public fluid_dynamics::InteractionWithWall diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp index 1749f6d746..4d91fae092 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp @@ -71,10 +71,6 @@ void BaseShearStressRelaxation2ndHalfWithWallvelocity_gradient_[index_i] += velocity_gradient; } -//============================================Plasticity===========================================// -//=================================================================================================// -//============BaseStressRelaxation1stHalfWithWall & BaseStressRelaxation2ndHalfWithWall============// -//=================================================================================================// template Vecd BaseStressRelaxation1stHalfWithWall::computeNonConservativeAcceleration(size_t index_i) { @@ -102,12 +98,10 @@ void BaseStressRelaxation1stHalfWithWall::inter Vecd& e_ij = wall_neighborhood.e_ij_[n]; Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n]; Real r_ij = wall_neighborhood.r_ij_[n]; - Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(Real(0), face_wall_external_acceleration); - Matd stress_tensor_in_wall = stress_tensor_i; - acceleration += (stress_tensor_i + stress_tensor_in_wall) * nablaW_ijV_j; + acceleration += 2 * stress_tensor_i * wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; } } @@ -138,10 +132,7 @@ void BaseStressRelaxation2ndHalfWithWall::inter Vecd& e_ij = wall_neighborhood.e_ij_[n]; Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n]; Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n]; - - Real beta = 1.9; - Vecd vel_in_wall = (1 - beta) * vel_i + beta * vel_ave_k[index_j]; - + Vecd vel_in_wall = 1.9 * vel_ave_k[index_j] - 0.9 *vel_i; Matd velocity_gradient_ij = -(vel_i - vel_in_wall) * nablaW_ijV_j.transpose(); velocity_gradient += velocity_gradient_ij; @@ -154,9 +145,7 @@ void BaseStressRelaxation2ndHalfWithWall::inter this->velocity_gradient_[index_i] += velocity_gradient; this->acc_[index_i] += p_dissipation / this->rho_[index_i]; } -//=================================================================================================// - //================================BaseStressDiffusionWithWall======================================// - //=================================================================================================// + template void BaseStressDiffusionWithWall::interaction(size_t index_i, Real dt) {