diff --git a/CMakeLists.txt b/CMakeLists.txt index 705feb1586..7809420646 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,10 +54,11 @@ target_link_libraries(sphinxsys_core INTERFACE ${Simbody_LIBS}) target_compile_definitions(sphinxsys_core INTERFACE _SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING) # Because of Simbody headers ### Eigen find_package(Eigen3 CONFIG REQUIRED) -target_link_libraries(sphinxsys_core INTERFACE Eigen3::Eigen) +target_link_libraries (sphinxsys_core INTERFACE Eigen3::Eigen) ### TBB find_package(TBB CONFIG REQUIRED) target_link_libraries(sphinxsys_core INTERFACE TBB::tbb TBB::tbbmalloc $<$:TBB::tbbmalloc_proxy>) +add_compile_definitions(TBB_SUPPRESS_DEPRECATED_MESSAGES=1) ### Threads find_package(Threads REQUIRED) target_link_libraries(sphinxsys_core INTERFACE Threads::Threads) diff --git a/PythonScriptStore/RegressionTest/regression_test_base_tool.py b/PythonScriptStore/RegressionTest/regression_test_base_tool.py index 1a997c6282..2fa1132eb8 100644 --- a/PythonScriptStore/RegressionTest/regression_test_base_tool.py +++ b/PythonScriptStore/RegressionTest/regression_test_base_tool.py @@ -16,22 +16,22 @@ def __init__(self, casename, bodyname, parametername): self.sphinxsys_case_name = casename self.sphinxsys_body_name = bodyname self.sphinxsys_parameter_name = parametername - self.enter_sphinxsys_exec_folder = f"cd {self.sphinxsys_exec_path};" - self.enter_sphinxsys_case_folder = f"cd {self.sphinxsys_case_path};" + self.enter_sphinxsys_exec_folder = f"{self.sphinxsys_exec_path}" + self.enter_sphinxsys_case_folder = f"{self.sphinxsys_case_path}" self.input_file_path = os.path.join(self.sphinxsys_exec_path, "input") self.condition_file_path = os.path.join(self.input_file_path, f"{bodyname}_{parametername}_runtimes.dat") def compile_case(self) -> None: print('Start compiling test case....') command = "make -j8" - os.system(self.enter_sphinxsys_case_folder) + os.chdir(self.enter_sphinxsys_case_folder) os.system(command) print('Compiling test case is finished...') def run_particle_relaxation(self) -> None: print('Start particle relaxation for the simulation...') command = f".{os.sep}{self.sphinxsys_case_name} --r=true" - os.system(self.enter_sphinxsys_exec_folder) + os.chdir(self.enter_sphinxsys_exec_folder) os.system(command) print('Simulating case is finished...') @@ -39,7 +39,7 @@ def run_case(self) -> None: print('Start case simulation...') print(self.enter_sphinxsys_exec_folder) command = f".{os.sep}{self.sphinxsys_case_name} --r=false --i=true --rt=true" - os.system(self.enter_sphinxsys_exec_folder) + os.chdir(self.enter_sphinxsys_exec_folder) os.system(command) print('Simulating case is finished...') @@ -60,36 +60,36 @@ def __init__(self, casename, bodyname, parametername): self.sphinxsys_case_name = casename self.sphinxsys_body_name = bodyname self.sphinxsys_parameter_name = parametername - self.enter_sphinxsys_exec_folder = f"cd {self.sphinxsys_exec_path};" - self.enter_sphinxsys_case_folder = f"cd {self.sphinxsys_case_path};" + self.enter_sphinxsys_exec_folder = f"{self.sphinxsys_exec_path}" + self.enter_sphinxsys_case_folder = f"{self.sphinxsys_case_path}" self.input_file_path = os.path.join(self.sphinxsys_exec_path, "input") self.condition_file_path = os.path.join(self.input_file_path, f"{bodyname}_{parametername}_runtimes.dat") def compile_case(self) -> None: print('Start compiling test case....') command = "make -j8" - os.system(self.enter_sphinxsys_case_folder) + os.chdir(self.enter_sphinxsys_case_folder) os.system(command) print('Compiling test case is finished...') def test_case(self) -> None: print('Start test case...') command = "make test" - os.system(self.enter_sphinxsys_case_folder) + os.chdir(self.enter_sphinxsys_case_folder) os.system(command) print('Testing case is finished...') def copy_reload(self) -> None: print('Start copy the reload file...') command = "cp -r reload bin" - os.system(self.enter_sphinxsys_case_folder) + os.chdir(self.enter_sphinxsys_case_folder) os.system(command) print('Copying the reload file is finished...') def run_particle_relaxation(self) -> None: print('Start particle relaxation for the simulation...') command = f".{os.sep}{self.sphinxsys_case_name} --r=true" - os.system(self.enter_sphinxsys_exec_folder) + os.chdir(self.enter_sphinxsys_exec_folder) os.system(command) print('Simulating case is finished...') @@ -97,7 +97,7 @@ def run_case(self) -> None: print('Start case simulation...') print(self.enter_sphinxsys_exec_folder) command = f".{os.sep}{self.sphinxsys_case_name} --r=false --i=true --rt=true" - os.system(self.enter_sphinxsys_exec_folder) + os.chdir(self.enter_sphinxsys_exec_folder) os.system(command) print('Simulating case is finished...') diff --git a/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.h b/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.h index 01712abab8..5fc8ed8463 100644 --- a/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.h +++ b/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.h @@ -182,11 +182,36 @@ namespace SPH void interaction(size_t index_i, Real dt = 0.0); private: - Real eta_; /**< damping coefficient */ StdLargeVec &Vol_, &mass_; StdLargeVec &variable_; StdVec *> wall_Vol_; StdVec *> wall_variable_; + + protected: + Real eta_; /**< damping coefficient */ + }; + + /** + * @class DampingPairwiseFromShell + * @brief Damping to wall by which the wall velocity is not updated + * and the mass of wall particle is not considered. + */ + template + class DampingPairwiseFromShell : public LocalDynamics, + public DataDelegateContact + { + public: + DampingPairwiseFromShell(BaseContactRelation &contact_relation, const std::string &variable_name, Real eta); + virtual ~DampingPairwiseFromShell(){}; + void interaction(size_t index_i, Real dt = 0.0); + + private: + StdLargeVec &Vol_, &mass_; + StdLargeVec &variable_; + StdVec *> shell_variable_; + + protected: + Real eta_; /**< damping coefficient */ }; /** diff --git a/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp b/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp index b83deeab59..cc1ccccbfc 100644 --- a/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp +++ b/SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp @@ -427,6 +427,53 @@ namespace SPH } } //=================================================================================================// + template + DampingPairwiseFromShell:: + DampingPairwiseFromShell(BaseContactRelation &contact_relation, const std::string &variable_name, Real eta) + : LocalDynamics(contact_relation.sph_body_), + DataDelegateContact(contact_relation), + eta_(eta), Vol_(particles_->Vol_), mass_(particles_->mass_), + variable_(*particles_->getVariableByName(variable_name)) + { + shell_variable_.reserve(contact_particles_.size()); + for (const auto& cp: contact_particles_) shell_variable_.push_back(cp->template getVariableByName(variable_name)); + } + //=================================================================================================// + template + void DampingPairwiseFromShell::interaction(size_t index_i, Real dt) + { + Real Vol_i = Vol_[index_i]; + Real mass_i = mass_[index_i]; + VariableType &variable_i = variable_[index_i]; + + std::array parameter_b; + + /** Contact interaction. */ + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + StdLargeVec &variable_k = *(shell_variable_[k]); + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + // forward sweep + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + + parameter_b[n] = contact_particles_[k]->DegeneratedSpacing(index_j) * eta_ * contact_neighborhood.dW_ijV_j_[n] * Vol_i * dt / contact_neighborhood.r_ij_[n]; + + // only update particle i + variable_[index_i] += parameter_b[n] * (variable_i - variable_k[index_j]) / (mass_i - 2.0 * parameter_b[n]); + } + // backward sweep + for (size_t n = contact_neighborhood.current_size_; n != 0; --n) + { + size_t index_j = contact_neighborhood.j_[n - 1]; + + // only update particle i + variable_[index_i] += parameter_b[n - 1] * (variable_i - variable_k[index_j]) / (mass_i - 2.0 * parameter_b[n - 1]); + } + } + } + //=================================================================================================// template template DampingWithRandomChoice:: diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h index a51314148b..00589a4a4a 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h @@ -40,3 +40,5 @@ #include "fluid_boundary.h" #include "all_eulerian_compressible_fluid_dynamics.h" #include "all_eulerian_weakly_compressible_fluid_dynamics.h" +#include "fluid_shell_interaction.h" +#include "fluid_shell_interaction.hpp" diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.cpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.cpp index 04e94f568b..8dd92cf34a 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.cpp +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.cpp @@ -36,7 +36,7 @@ namespace SPH Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; + acceleration_trans -= 2.0 * nablaW_ijV_j * this->contact_particles_[k]->DegeneratedSpacing(index_j); } } @@ -59,7 +59,7 @@ namespace SPH Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n]; // acceleration for transport velocity - acceleration_trans -= 2.0 * nablaW_ijV_j; + acceleration_trans -= 2.0 * nablaW_ijV_j * this->contact_particles_[k]->DegeneratedSpacing(index_j); } } diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp index 87a28bb644..b520f23ba8 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp @@ -88,7 +88,8 @@ namespace SPH Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { - sigma += contact_neighborhood.W_ij_[n] * contact_inv_rho0_k * contact_mass_k[contact_neighborhood.j_[n]]; + Real mass = this->contact_particles_[k]->ParticleMass(contact_neighborhood.j_[n]); + sigma += contact_neighborhood.W_ij_[n] * contact_inv_rho0_k * mass; } } return sigma; diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp index 90375c5fa0..adf5311534 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp @@ -14,6 +14,7 @@ namespace SPH //=================================================================================================// BaseDensitySummationInner::BaseDensitySummationInner(BaseInnerRelation &inner_relation) : LocalDynamics(inner_relation.sph_body_), FluidDataInner(inner_relation), + spacing_ref_(sph_body_.sph_adaptation_->ReferenceSpacing()), rho_(particles_->rho_), rho_sum_(particles_->rho_sum_), mass_(particles_->mass_), rho0_(sph_body_.base_material_->ReferenceDensity()) {} //=================================================================================================// diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h index 0f05d60312..156d82dcec 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h @@ -75,6 +75,7 @@ namespace SPH protected: Real rho0_; + Real spacing_ref_; StdLargeVec &rho_, &rho_sum_, &mass_; }; diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.cpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.cpp new file mode 100644 index 0000000000..9b7d9b10bc --- /dev/null +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.cpp @@ -0,0 +1,12 @@ +#include "fluid_shell_interaction.h" + +//=================================================================================================// +namespace SPH +{ + //=================================================================================================// + namespace fluid_dynamics + { + //=================================================================================================// + } +//=================================================================================================// +} \ No newline at end of file diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.h b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.h new file mode 100644 index 0000000000..80862df59c --- /dev/null +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.h @@ -0,0 +1,136 @@ +/* -------------------------------------------------------------------------* + * 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-2022 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_shell_interaction.h +* @brief Here, we define the algorithm classes for the interaction between fluid and +* thin structure, plate and shell. +* @author Chi ZHang and Xiangyu Hu +*/ + + +#ifndef FLUID_SHELL_INTERACTION_H +#define FLUID_SHELL_INTERACTION_H + +#include "fluid_dynamics_complex.h" +#include "fluid_dynamics_complex.hpp" + +namespace SPH +{ + namespace fluid_dynamics + { + typedef DataDelegateContact FluidShellContactData; + typedef DataDelegateContact FluidShellData; + + /** + * @class InteractionWithShell + * @brief Abstract base class for general fluid-shell interaction model. + * @note Here, + */ + template + class InteractionWithShell : public BaseIntegrationType, public FluidShellData + { + public: + template + InteractionWithShell(BaseContactRelation &contact_relation, BaseBodyRelationType &base_body_relation, Args &&...args); + + template + InteractionWithShell(ComplexRelation &fluid_shell_relation, Args &&...args) + : InteractionWithShell(fluid_shell_relation.getContactRelation(), fluid_shell_relation.getInnerRelation(), std::forward(args)...) + {} + + virtual ~InteractionWithShell(){}; + + protected: + Real spacing_ref_; + StdVec *> shell_thickness_; + StdVec *> shell_n_, shell_vel_ave_, shell_acc_ave_; + }; + + /** + * @class BaseFluidShellIntegration1stHalf + * @brief template class for fluid-shell pressure relaxation scheme + */ + template + class BaseFluidShellIntegration1stHalf : public InteractionWithShell + { + public: + template + explicit BaseFluidShellIntegration1stHalf(Args &&...args) + : InteractionWithShell(std::forward(args)...){}; + virtual ~BaseFluidShellIntegration1stHalf(){}; + void interaction(size_t index_i, Real dt = 0.0); + + protected: + virtual Vecd computeNonConservativeAcceleration(size_t index_i) override; + }; + + using FluidShellIntegration1stHalf = BaseFluidShellIntegration1stHalf; + using FluidShellIntegration1stHalfRiemann = BaseFluidShellIntegration1stHalf; + + using FluidShellandWallIntegration1stHalf = BaseIntegration1stHalfWithWall; + using FluidShellandWallIntegration1stHalfRiemann = BaseIntegration1stHalfWithWall; + using ExtendFluidShellandWallIntegration1stHalfRiemann = BaseExtendIntegration1stHalfWithWall; + + /** + * @class BaseFluidShellIntegration2ndHalf + * @brief template class pressure relaxation scheme with wall boundary + */ + template + class BaseFluidShellIntegration2ndHalf : public InteractionWithShell + { + public: + template + explicit BaseFluidShellIntegration2ndHalf(Args &&...args) + : InteractionWithShell(std::forward(args)...) + {}; + virtual ~BaseFluidShellIntegration2ndHalf(){}; + void interaction(size_t index_i, Real dt = 0.0); + }; + using FluidShellIntegration2ndHalf = BaseFluidShellIntegration2ndHalf; + using FluidShellIntegration2ndHalfRiemann = BaseFluidShellIntegration2ndHalf; + + using FluidShellandWallIntegration2ndHalf = BaseIntegration2ndHalfWithWall; + using FluidShellandWallIntegration2ndHalfRiemann = BaseIntegration2ndHalfWithWall; + + /** + * @class ViscousWithWall + * @brief template class viscous acceleration with wall boundary + */ + template + class BaseViscousAccelerationWithShell : public InteractionWithShell + { + public: + template + BaseViscousAccelerationWithShell(Args &&...args) + : InteractionWithShell(std::forward(args)...) + {}; + virtual ~BaseViscousAccelerationWithShell(){}; + void interaction(size_t index_i, Real dt = 0.0); + }; + + using ViscousAccelerationWithShell = BaseViscousAccelerationWithShell; + using ViscousAccelerationWithShellandWall = BaseViscousAccelerationWithWall; + + } +} +#endif //FLUID_SHELL_INTERACTION_H \ No newline at end of file diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.hpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.hpp new file mode 100644 index 0000000000..fbe86411ce --- /dev/null +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.hpp @@ -0,0 +1,171 @@ +/* -------------------------------------------------------------------------* + * 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-2022 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_shell_interaction.h +* @brief Here, we define the algorithm classes for the interaction between fluid and +* thin structure, plate and shell. +* @author Chi ZHang and Xiangyu Hu +*/ +#pragma once + +#include "fluid_shell_interaction.h" + +//=================================================================================================// +namespace SPH +{ +//=================================================================================================// + namespace fluid_dynamics + { + //=================================================================================================// + template + template + InteractionWithShell::InteractionWithShell(BaseContactRelation &contact_relation, BaseBodyRelationType &base_body_relation, Args &&...args) + : BaseIntegrationType(base_body_relation, std::forward(args)...) + , FluidShellData(contact_relation) + { + if (&base_body_relation.sph_body_ != &contact_relation.sph_body_) + { + std::cout << "\n Error: the two body_realtions do not have the same source body!" << std::endl; + std::cout << __FILE__ << ':' << __LINE__ << std::endl; + exit(1); + } + + shell_thickness_.reserve(FluidShellData::contact_particles_.size()); + shell_n_.reserve(FluidShellData::contact_particles_.size()); + shell_vel_ave_.reserve(FluidShellData::contact_particles_.size()); + shell_acc_ave_.reserve(FluidShellData::contact_particles_.size()); + + for (const auto& cp: FluidShellData::contact_particles_) + { + shell_thickness_.push_back(&(cp->thickness_)); + shell_n_.push_back(&(cp->n_)); + shell_vel_ave_.push_back(cp->AverageVelocity()); + shell_acc_ave_.push_back(cp->AverageAcceleration()); + } + } + //=================================================================================================// + template + void BaseFluidShellIntegration1stHalf::interaction(size_t index_i, Real dt) + { + Integration1stHalfType::interaction(index_i, dt); + Vecd acc_prior_i = computeNonConservativeAcceleration(index_i); + + Vecd acceleration = Vecd::Zero(); + Real rho_dissipation = 0.0; + for (size_t k = 0; k < FluidShellData::contact_configuration_.size(); ++k) + { + StdLargeVec &vel_ave_k = *(this->shell_vel_ave_[k]); + StdLargeVec &acc_ave_k = *(this->shell_acc_ave_[k]); + StdLargeVec &n_k = *(this->shell_n_[k]); + StdLargeVec &thickness_k = *(this->shell_thickness_[k]); + + Neighborhood &shell_neighborhood = (*FluidShellData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != shell_neighborhood.current_size_; ++n) + { + size_t index_j = shell_neighborhood.j_[n]; + Vecd &e_ij = shell_neighborhood.e_ij_[n]; + Real dW_ijV_j = shell_neighborhood.dW_ijV_j_[n]; + Real r_ij = shell_neighborhood.r_ij_[n]; + + Real face_shell_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); + Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration); + acceleration -= (this->p_[index_i] + p_in_shell) * e_ij * dW_ijV_j * thickness_k[index_j]; + rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_shell) * dW_ijV_j * thickness_k[index_j]; + } + } + this->acc_[index_i] += acceleration / this->rho_[index_i]; + this->drho_dt_[index_i] += rho_dissipation * this->rho_[index_i]; + } + //=================================================================================================// + template + Vecd BaseFluidShellIntegration1stHalf::computeNonConservativeAcceleration(size_t index_i) + { + return this->acc_prior_[index_i]; + } + //=================================================================================================// + template + void BaseFluidShellIntegration2ndHalf::interaction(size_t index_i, Real dt) + { + BaseIntegration2ndHalfType::interaction(index_i, dt); + + Real density_change_rate = 0.0; + Vecd p_dissipation = Vecd::Zero(); + + for (size_t k = 0; k < FluidShellData::contact_configuration_.size(); ++k) + { + Vecd &acc_prior_i = this->acc_prior_[index_i]; + + StdLargeVec &vel_ave_k = *(this->shell_vel_ave_[k]); + StdLargeVec &acc_ave_k = *(this->shell_acc_ave_[k]); + StdLargeVec &n_k = *(this->shell_n_[k]); + StdLargeVec &thickness_k = *(this->shell_thickness_[k]); + Neighborhood &shell_neighborhood = (*FluidShellData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != shell_neighborhood.current_size_; ++n) + { + size_t index_j = shell_neighborhood.j_[n]; + Vecd &e_ij = shell_neighborhood.e_ij_[n]; + Real r_ij = shell_neighborhood.r_ij_[n]; + Real dW_ijV_j = shell_neighborhood.dW_ijV_j_[n]; + Vecd correct_n = SGN( e_ij.dot(n_k[index_j]) ) * n_k[index_j]; + + Vecd vel_in_shell = 2.0 * vel_ave_k[index_j] - this->vel_[index_i]; + Real u_jump = (this->vel_[index_i] - vel_in_shell).dot(correct_n); + density_change_rate += u_jump * dW_ijV_j * thickness_k[index_j]; + p_dissipation += this->riemann_solver_.DissipativePJump(u_jump) * dW_ijV_j * thickness_k[index_j] * e_ij; + } + } + this->drho_dt_[index_i] += density_change_rate * this->rho_[index_i]; + this->acc_[index_i] += p_dissipation / this->rho_[index_i]; + } + //=================================================================================================// + template + void BaseViscousAccelerationWithShell::interaction(size_t index_i, Real dt) + { + ViscousAccelerationInnerType::interaction(index_i, dt); + + Real rho_i = this->rho_[index_i]; + const Vecd &vel_i = this->vel_[index_i]; + + Vecd acceleration = Vecd::Zero(); + Vecd vel_derivative = Vecd::Zero(); + for (size_t k = 0; k < FluidShellData::contact_configuration_.size(); ++k) + { + StdLargeVec &vel_ave_k = *(this->shell_vel_ave_[k]); + StdLargeVec &thickness_k = *(this->shell_thickness_[k]); + Neighborhood &contact_neighborhood = (*FluidShellData::contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Real r_ij = contact_neighborhood.r_ij_[n]; + + vel_derivative = 2.0 * (vel_i - vel_ave_k[index_j]) / (r_ij + 0.01 * this->smoothing_length_); + acceleration += 2.0 * this->mu_ * vel_derivative * contact_neighborhood.dW_ijV_j_[n] * thickness_k[index_j] / rho_i; + } + } + + this->acc_prior_[index_i] += acceleration; + } + //=================================================================================================// + } +//=================================================================================================// +} \ No newline at end of file diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp index 03047b895e..6ac11ed87c 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.cpp @@ -6,24 +6,6 @@ namespace SPH namespace fluid_dynamics { //=================================================================================================// - FreeSurfaceIndicationComplex:: - FreeSurfaceIndicationComplex(BaseInnerRelation &inner_relation, - BaseContactRelation &contact_relation, Real threshold) - : FreeSurfaceIndicationInner(inner_relation, threshold), FluidContactData(contact_relation) - { - for (size_t k = 0; k != contact_particles_.size(); ++k) - { - Real rho0_k = contact_bodies_[k]->base_material_->ReferenceDensity(); - contact_inv_rho0_.push_back(1.0 / rho0_k); - contact_mass_.push_back(&(contact_particles_[k]->mass_)); - } - } - //=================================================================================================// - FreeSurfaceIndicationComplex:: - FreeSurfaceIndicationComplex(ComplexRelation &complex_relation, Real threshold) - : FreeSurfaceIndicationComplex(complex_relation.getInnerRelation(), - complex_relation.getContactRelation(), threshold) {} - //=================================================================================================// void FreeSurfaceIndicationComplex::interaction(size_t index_i, Real dt) { FreeSurfaceIndicationInner::interaction(index_i, dt); @@ -34,7 +16,8 @@ namespace SPH 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 -= contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.r_ij_[n] + * this->contact_particles_[k]->DegeneratedSpacing(contact_neighborhood.j_[n]); } } pos_div_[index_i] += pos_div; diff --git a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h index c78f132c22..b03ba6e6b9 100644 --- a/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h +++ b/SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_surface_complex.h @@ -39,24 +39,21 @@ namespace SPH namespace fluid_dynamics { /** - * @class FreeSurfaceIndicationComplex - * @brief indicate the particles near the free fluid surface. + * @class FreeSurfaceIndicationComplex + * @brief indicate the particles near the free fluid surface. */ - class FreeSurfaceIndicationComplex : public FreeSurfaceIndicationInner, public FluidContactData + class FreeSurfaceIndicationComplex + : public BaseInteractionComplex { public: - FreeSurfaceIndicationComplex(BaseInnerRelation &inner_relation, - BaseContactRelation &contact_relation, Real threshold = 0.75); - explicit FreeSurfaceIndicationComplex(ComplexRelation &complex_relation, Real threshold = 0.75); + template + FreeSurfaceIndicationComplex(Args &&...args) + : BaseInteractionComplex(std::forward(args)...){}; virtual ~FreeSurfaceIndicationComplex(){}; void interaction(size_t index_i, Real dt = 0.0); - - protected: - StdVec contact_inv_rho0_; - StdVec *> contact_mass_; }; - using SpatialTemporalFreeSurfaceIdentificationComplex = - SpatialTemporalFreeSurfaceIdentification; + + using SpatialTemporalFreeSurfaceIdentificationComplex =SpatialTemporalFreeSurfaceIdentification; /** the cases with free surface and freestream */ using DensitySummationFreeSurfaceComplex = DensitySummationFreeSurface; diff --git a/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h index 123d99212a..b9dd1fe91e 100644 --- a/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h +++ b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h @@ -38,3 +38,4 @@ #include "fluid_structure_interaction.h" #include "thin_structure_math.h" #include "inelastic_dynamics.h" +#include "shell_fluid_interaction.h" diff --git a/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.cpp b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.cpp new file mode 100644 index 0000000000..e1ab8739d7 --- /dev/null +++ b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.cpp @@ -0,0 +1,86 @@ +#include "shell_fluid_interaction.h" + +namespace SPH +{ + //=====================================================================================================// + namespace solid_dynamics + { + //=================================================================================================// + FluidViscousForceOnShell::FluidViscousForceOnShell(BaseContactRelation &contact_relation) + : LocalDynamics(contact_relation.sph_body_) + , FluidShellContactData(contact_relation) + , Vol_(particles_->Vol_) + , vel_ave_(*particles_->AverageVelocity()) + { + particles_->registerVariable(viscous_force_from_fluid_, "ViscousForceFromFluid"); + + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_fluids_.push_back(&contact_particles_[k]->fluid_); + mu_.push_back(contact_fluids_[k]->ReferenceViscosity()); + smoothing_length_.push_back(contact_bodies_[k]->sph_adaptation_->ReferenceSmoothingLength()); + + contact_rho_n_.push_back(&(contact_particles_[k]->rho_)); + contact_vel_n_.push_back(&(contact_particles_[k]->vel_)); + } + } + //=================================================================================================// + void FluidViscousForceOnShell::interaction(size_t index_i, Real dt) + { + Real Vol_i = Vol_[index_i]; + const Vecd &vel_ave_i = vel_ave_[index_i]; + + Vecd force = Vecd::Zero(); + /** Contact interaction. */ + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Real mu_k = mu_[k]; + Real smoothing_length_k = smoothing_length_[k]; + StdLargeVec &vel_n_k = *(contact_vel_n_[k]); + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + + Vecd vel_derivative = 2.0 * (vel_ave_i - vel_n_k[index_j]) / + (contact_neighborhood.r_ij_[n] + 0.01 * smoothing_length_k); + + force += 2.0 * mu_k * vel_derivative * contact_neighborhood.dW_ijV_j_[n] * Vol_i * particles_->DegeneratedSpacing(index_i); + } + } + + viscous_force_from_fluid_[index_i] = force; + } + //=================================================================================================// + void FluidAngularConservativeViscousForceOnShell::interaction(size_t index_i, Real dt) + { + Real Vol_i = Vol_[index_i]; + const Vecd &vel_ave_i = vel_ave_[index_i]; + + Vecd force = Vecd::Zero(); + /** Contact interaction. */ + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Real mu_k = mu_[k]; + Real smoothing_length_k = smoothing_length_[k]; + StdLargeVec &rho_n_k = *(contact_rho_n_[k]); + StdLargeVec &vel_n_k = *(contact_vel_n_[k]); + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + + /** The following viscous force is given in Monaghan 2005 (Rep. Prog. Phys.) */ + Real v_r_ij = contact_neighborhood.r_ij_[n] * (vel_ave_i - vel_n_k[index_j]).dot(contact_neighborhood.e_ij_[n]); + Real vel_difference = 0.0 * (vel_ave_i - vel_n_k[index_j]).norm() * contact_neighborhood.r_ij_[n]; + Real eta_ij = 8.0 * SMAX(mu_k, rho_n_k[index_j] * vel_difference) * v_r_ij / + (contact_neighborhood.r_ij_[n] * contact_neighborhood.r_ij_[n] + 0.01 * smoothing_length_k); + force += eta_ij * contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n] * Vol_i * particles_->DegeneratedSpacing(index_i); + } + } + + viscous_force_from_fluid_[index_i] = force; + } + //=================================================================================================// + } +} diff --git a/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.h b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.h new file mode 100644 index 0000000000..74fbe80678 --- /dev/null +++ b/SPHINXsys/src/shared/particle_dynamics/solid_dynamics/shell_fluid_interaction.h @@ -0,0 +1,207 @@ +/* -------------------------------------------------------------------------* + * 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-2022 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 shell_fluid_interaction.h + * @brief Here, we define the algorithm classes for fluid force on shell structure. + * @author Chi Zhang and Xiangyu Hu + */ + +#ifndef SHELL_FLUID_INTERACTION_H +#define SHELL_FLUID_INTERACTION_H + +#include "all_particle_dynamics.h" +#include "base_material.h" +#include "fluid_dynamics_complex.h" +#include "elastic_dynamics.h" +#include "riemann_solver.h" + +namespace SPH +{ + namespace solid_dynamics + { + typedef DataDelegateContact FluidShellContactData; + /** + * @class FluidViscousForceOnShell + * @brief Computing the viscous force on shell structure from the fluid + */ + class FluidViscousForceOnShell : public LocalDynamics, public FluidShellContactData + { + public: + explicit FluidViscousForceOnShell(BaseContactRelation &contact_relation); + virtual ~FluidViscousForceOnShell(){}; + + void interaction(size_t index_i, Real dt = 0.0); + StdLargeVec &getViscousForceFromFluid() { return viscous_force_from_fluid_; }; + + protected: + StdLargeVec &Vol_; + StdLargeVec &vel_ave_; + StdLargeVec viscous_force_from_fluid_; + + StdVec contact_fluids_; + StdVec mu_; + StdVec smoothing_length_; + + StdVec *> contact_rho_n_; + StdVec *> contact_vel_n_; + }; + + /** + * @class FluidAngularConservativeViscousForceOnShell + * @brief Computing the viscous force on shell structure from the fluid + * TODO: new test for this. + */ + class FluidAngularConservativeViscousForceOnShell : public FluidViscousForceOnShell + { + public: + explicit FluidAngularConservativeViscousForceOnShell(BaseContactRelation &contact_relation) + : FluidViscousForceOnShell(contact_relation) + {}; + virtual ~FluidAngularConservativeViscousForceOnShell(){}; + + protected: + void interaction(size_t index_i, Real dt = 0.0); + }; + + /** + * @class BaseFluidPressureForceOnShell + * @brief Template class fro computing the pressure force from the fluid with different Riemann solvers. + * The pressure force is added on the viscous force of the latter is computed. + * This class is for fluid-shell-interaction applications to achieve smaller solid dynamics + * time step size compared to the fluid dynamics + */ + template + class BaseFluidPressureForceOnShell : public LocalDynamics, public FluidShellContactData + { + public: + explicit BaseFluidPressureForceOnShell(BaseContactRelation &contact_relation) + : LocalDynamics(contact_relation.sph_body_) + , FluidShellContactData(contact_relation) + , vel_ave_(*particles_->AverageVelocity()) + , acc_prior_(particles_->acc_prior_) + , acc_ave_(*particles_->AverageAcceleration()) + , n_(particles_->n_) + { + particles_->registerVariable(force_from_fluid_, "ForceFromFluid"); + + contact_fluids_.reserve(contact_particles_.size()); + riemann_solvers_.reserve(contact_particles_.size()); + contact_rho_n_.reserve(contact_particles_.size()); + contact_p_.reserve(contact_particles_.size()); + contact_vel_n_.reserve(contact_particles_.size()); + contact_acc_prior_.reserve(contact_particles_.size()); + + for (const auto& cp: contact_particles_) + { + contact_fluids_.push_back(&cp->fluid_); + riemann_solvers_.push_back(RiemannSolverType(cp->fluid_, cp->fluid_)); + + contact_rho_n_.push_back(&(cp->rho_)); + contact_p_.push_back(&(cp->p_)); + contact_vel_n_.push_back(&(cp->vel_)); + contact_acc_prior_.push_back(&(cp->acc_prior_)); + } + + }; + virtual ~BaseFluidPressureForceOnShell(){}; + + void interaction(size_t index_i, Real dt = 0.0) + { + const Vecd &acc_ave_i = acc_ave_[index_i]; + const Vecd &vel_ave_i = vel_ave_[index_i]; + const Vecd &n_i = n_[index_i]; + + Vecd force = Vecd::Zero(); + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + StdLargeVec &rho_n_k = *(contact_rho_n_[k]); + StdLargeVec &p_k = *(contact_p_[k]); + StdLargeVec &vel_n_k = *(contact_vel_n_[k]); + StdLargeVec &acc_prior_k = *(contact_acc_prior_[k]); + Fluid *fluid_k = contact_fluids_[k]; + RiemannSolverType &riemann_solver_k = riemann_solvers_[k]; + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Vecd e_ij = contact_neighborhood.e_ij_[n]; + Real r_ij = contact_neighborhood.r_ij_[n]; + Vecd correct_n = -SGN( e_ij.dot(n_i) ) * n_i; + + Real face_wall_external_acceleration = (acc_prior_k[index_j] - acc_ave_i).dot(e_ij); + Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(0.0, face_wall_external_acceleration); + + Vecd vel_in_shell = 2.0 * vel_ave_i - vel_n_k[index_j]; + Real u_jump = (vel_n_k[index_j] - vel_in_shell).dot(correct_n); + + force -= (p_in_wall + p_k[index_j] - riemann_solver_k.DissipativePJump(u_jump)) * e_ij + * contact_neighborhood.dW_ijV_j_[n] * particles_->ParticleVolume(index_i); + } + } + force_from_fluid_[index_i] = force; + acc_prior_[index_i] = force / particles_->ParticleMass(index_i); + }; + + protected: + StdLargeVec &vel_ave_, &acc_prior_, &acc_ave_, &n_; + StdLargeVec force_from_fluid_; + + StdVec contact_fluids_; + StdVec riemann_solvers_; + StdVec *> contact_rho_n_, contact_p_; + StdVec *> contact_vel_n_, contact_acc_prior_; + }; + using FluidPressureForceOnShell = BaseFluidPressureForceOnShell; + using FluidPressureForceOnShellRiemann = BaseFluidPressureForceOnShell; + + /** + * @class BaseFluidForceOnShellUpdate + * @brief template class for computing force on shell structure from fluid with updated viscous force + */ + template + class BaseFluidForceOnShellUpdate : public PressureForceType + { + public: + template + BaseFluidForceOnShellUpdate(BaseContactRelation &contact_relation, ViscousForceOnShellType &viscous_force_on_solid) + : PressureForceType(contact_relation) + , viscous_force_from_fluid_(viscous_force_on_solid.getViscousForceFromFluid()) + {}; + virtual ~BaseFluidForceOnShellUpdate(){}; + + void interaction(size_t index_i, Real dt = 0.0) + { + PressureForceType::interaction(index_i, dt); + this->force_from_fluid_[index_i] += viscous_force_from_fluid_[index_i]; + this->acc_prior_[index_i] += viscous_force_from_fluid_[index_i] / this->particles_->ParticleMass(index_i); + }; + + protected: + StdLargeVec &viscous_force_from_fluid_; + }; + + using FluidForceOnShellUpdate = BaseFluidForceOnShellUpdate; + using FluidForceOnShellUpdateRiemann = BaseFluidForceOnShellUpdate; + } +} +#endif // FLUID_STRUCTURE_INTERACTION_H \ No newline at end of file diff --git a/SPHINXsys/src/shared/particles/base_particles.h b/SPHINXsys/src/shared/particles/base_particles.h index 2a18dc3f3c..4f0ae7873a 100644 --- a/SPHINXsys/src/shared/particles/base_particles.h +++ b/SPHINXsys/src/shared/particles/base_particles.h @@ -200,7 +200,8 @@ namespace SPH virtual Real ParticleVolume(size_t index_i) { return Vol_[index_i]; } /** Return particle mass. */ virtual Real ParticleMass(size_t index_i) { return mass_[index_i]; } - + /** Return degenerated spacing. */ + virtual Real DegeneratedSpacing(size_t index_i) { return 1.0; } protected: SPHBody &sph_body_; /**< The body in which the particles belongs to. */ std::string body_name_; /**< Name of the body. */ diff --git a/SPHINXsys/src/shared/particles/fluid_particles.cpp b/SPHINXsys/src/shared/particles/fluid_particles.cpp index 1105bdab61..1e75c841e4 100644 --- a/SPHINXsys/src/shared/particles/fluid_particles.cpp +++ b/SPHINXsys/src/shared/particles/fluid_particles.cpp @@ -32,6 +32,7 @@ namespace SPH // add restart output particle data //---------------------------------------------------------------------- addVariableToRestart("Pressure"); + addVariableToWrite("Density"); } //=================================================================================================// ViscoelasticFluidParticles:: diff --git a/SPHINXsys/src/shared/particles/solid_particles.h b/SPHINXsys/src/shared/particles/solid_particles.h index b7fe3c99d2..b4bfbaeba5 100644 --- a/SPHINXsys/src/shared/particles/solid_particles.h +++ b/SPHINXsys/src/shared/particles/solid_particles.h @@ -183,6 +183,13 @@ namespace SPH virtual Real ParticleVolume(size_t index_i) override { return Vol_[index_i] * thickness_[index_i]; } /** get particle mass. */ virtual Real ParticleMass(size_t index_i) override { return mass_[index_i] * thickness_[index_i]; } + /** Return degenerated spacing. */ + virtual Real DegeneratedSpacing(size_t index_i) override { return thickness_[index_i]; } + /** Get wall average velocity when interacting with fluid. */ + virtual StdLargeVec *AverageVelocity() override { return &vel_ave_; }; + /** Get wall average acceleration when interacting with fluid. */ + virtual StdLargeVec *AverageAcceleration() override { return &acc_ave_; }; + /** Initialize variable for shell particles. */ virtual void initializeOtherVariables() override; /** Return this pointer. */ diff --git a/cmake/FindSIMD.cmake b/cmake/FindSIMD.cmake index 5ab2984155..34be3a6913 100644 --- a/cmake/FindSIMD.cmake +++ b/cmake/FindSIMD.cmake @@ -49,7 +49,6 @@ function (test_sse_availability) __m128i va = _mm_loadu_si128((__m128i*)a); __m128i vb = _mm_loadu_si128((__m128i*)b); __m128i vc = _mm_cmpgt_epi64(va, vb); - _mm_storeu_si128((__m128i*)c, vc); if (c[0] == -1LL && c[1] == 0LL) return 0; @@ -75,7 +74,6 @@ function (test_sse_availability) __m128i va = _mm_loadu_si128((__m128i*)a); __m128i vb = _mm_loadu_si128((__m128i*)b); __m128i vc = _mm_cmpeq_epi64(va, vb); - _mm_storeu_si128((__m128i*)c, vc); if (c[0] == 0LL && c[1] == -1LL) return 0; @@ -96,17 +94,14 @@ function (test_sse_availability) #else #include #endif - int main() { float a[4] = { 1.0f, 2.0f, 3.0f, 4.0f }; float b[4] = { 3.0f, 5.0f, 7.0f, 9.0f }; float c[4]; - __m128 va = _mm_loadu_ps(a); __m128 vb = _mm_loadu_ps(b); __m128 vc = _mm_hadd_ps(va, vb); - _mm_storeu_ps(c, vc); if (c[0] == 3.0f && c[1] == 7.0f && c[2] == 8.0f && c[3] == 16.0f) return 0; @@ -209,7 +204,6 @@ function (test_fma_availability) __m256d a = _mm256_set_pd (-1, 2, -3, 4); __m256d b = _mm256_set_pd (-2, 3, -4, 1); __m256d c = _mm256_set_pd (-11, 6, 4, -1); - __m256d result = _mm256_fmsub_pd (a, b, c); return 0; }" DETECTED_FMA) @@ -377,12 +371,10 @@ function(test_neon_availability) float64_t a[2] = { 1., 2. }; float64_t b[2] = { -1., 3. }; float64_t c[2]; - float64x2_t va = vld1q_f64(&a[0]); float64x2_t vb = vld1q_f64(&b[0]); float64x2_t vc = vaddq_f64(va, vb); vst1q_f64(&c[0], vc); - if (c[0] == 0. && c[1] == 5.) return 0; else @@ -472,4 +464,4 @@ endif() set(SIMD_C_FLAGS "${SIMD_FLAGS}" CACHE STRING "Flags used for compiling C programs with SIMD support") set(SIMD_CXX_FLAGS "${SIMD_FLAGS}" CACHE STRING "Flags used for compiling C++ programs with SIMD support") -mark_as_advanced(SIMD_SSE SIMD_AVX SIMD_FMA SIMD_NEON SIMD_C_FLAGS SIMD_CXX_FLAGS) +mark_as_advanced(SIMD_SSE SIMD_AVX SIMD_FMA SIMD_NEON SIMD_C_FLAGS SIMD_CXX_FLAGS) \ No newline at end of file diff --git a/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp b/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp index b49846f9cc..3b658eca7c 100644 --- a/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp +++ b/tests/2d_examples/test_2d_eulerian_flow_around_cylinder/2d_eulerian_flow_around_cylinder.cpp @@ -168,7 +168,7 @@ int main(int ac, char *av[]) InteractionWithUpdate density_relaxation(water_block_complex); InteractionDynamics viscous_acceleration(water_block_complex); ReduceDynamics get_fluid_time_step_size(water_block); - InteractionWithUpdate surface_indicator(water_block_complex.getInnerRelation(), water_block_complex.getContactRelation()); + InteractionWithUpdate surface_indicator(water_block_complex); InteractionDynamics variable_reset_in_boundary_condition(water_block_complex.getInnerRelation()); //---------------------------------------------------------------------- // Compute the force exerted on solid body due to fluid pressure and viscosity diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/CMakeLists.txt b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/CMakeLists.txt new file mode 100644 index 0000000000..9c57fba9fd --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/CMakeLists.txt @@ -0,0 +1,23 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir + +set(CMAKE_VERBOSE_MAKEFILE on) + +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}) + +aux_source_directory(. DIR_SRCS) +add_executable(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS} ) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} sphinxsys_2d) + +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/case.h b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/case.h new file mode 100644 index 0000000000..3624c64d6f --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/case.h @@ -0,0 +1,227 @@ +/** +* @file case.h +* @brief Fluid-shell interaction in sloshing flow. +* @details Here, the first fluid-shell interaction test is presented. +* @author Chi Zhang +*/ + +#include "sphinxsys.h" +#define PI 3.1415926 + /** +* @brief Namespace cite here. +*/ +using namespace SPH; +/** + * @brief Basic geometry parameters and numerical setup. + */ +Real DL = 1.0; /**< Tank length. */ +Real DH = 0.7; /**< Tank height. */ +Real L_W = 1.0; /**< water width. */ +Real L_H = 0.15; /**< water depth. */ +Real Gate_x = 0.5 * L_W; /**< Width of the gate. */ +Real Gate_width = 0.006; /**< Width of the gate. */ +Real Gate_height = 0.174; /**< Height of the gate. */ +Real particle_spacing_ref = Gate_width; /**< Initial reference particle spacing. */ +Real BW = particle_spacing_ref * 4.0; /**< Extending width for BCs. */ +Real thickness = Gate_width; +int particle_number_mid_surface = (0.026 + Gate_height + BW) / particle_spacing_ref; +BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(DL + BW, DH + BW)); +// parameters for liquid sloshing Case:Xue Mian +Real a_0 = 0.01; /**< amplitude of the sloshing. */ +Real w_0 = 1.2 * 4.142814038; /**< frequency of the sloshing 0.583*8.9556. */ +Real k_0 = a_0 * w_0 * w_0; /**< parameter of sloshing x = k_0 * sin(2*pi*f_0*t). */ +/** + * @brief Material properties of the fluid. + */ +Real rho0_f = 1000.0; /**< Reference density of fluid. */ +Real gravity_g = 9.81; /**< Value of gravity. */ +Real U_max = 2.0*sqrt(gravity_g*L_H); /**< Characteristic velocity. */ +Real c_f = 10.0* U_max; /**< Reference sound speed. */ +Real mu_f = 1.0e-6; +/** + * @brief Material properties of the elastic gate. + */ +Real rho0_s = 1250; /**< Reference density of gate. */ +Real Youngs_modulus = 30.0e6; /**< reference Youngs modulus. */ +Real poisson = 0.47; /**< Poisson ratio. */ +Real physical_viscosity = 1000.0; /**< physical damping, here we choose the same value as numerical viscosity. */ +/** create a water block shape */ +StdVec water_block_shape{ + Vecd(0.0, 0.0), Vecd(0.0, L_H), Vecd(L_W, L_H), Vecd(L_W, 0.0), Vecd(0.0, 0.0)}; +/** Baffle base shape */ +StdVec baffle_base_shap{ + Vecd(Gate_x - 0.0262, 0.0), + Vecd(Gate_x - 0.0262, 0.006), + Vecd(Gate_x - 0.0262 + 0.0075, 0.006), + Vecd(Gate_x - 0.0262 + 0.0075, 0.0125), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008, 0.0125), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008, 0.026), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022, 0.026), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022, 0.0125), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022 + 0.008, 0.0125), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022 + 0.008, 0.006), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022 + 0.008 + 0.0075, 0.006), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.022 + 0.008 + 0.0075, 0.0), + Vecd(Gate_x - 0.0262, 0.0) +}; +/** Inner wall shape */ +StdVec inner_wall_shape{ + Vecd(0.0, 0.0), + Vecd(0.0, DH), + Vecd(DL, DH), + Vecd(DL, 0.0), + Vecd(0.0, 0.0) +}; +/** Outer wall shape */ +StdVec outer_wall_shape{ + Vecd(-BW, -BW), + Vecd(-BW, DH + BW), + Vecd(DL + BW, DH + BW), + Vecd(DL + BW, -BW), + Vecd(-BW, -BW) +}; +/** Baffle shape. */ +StdVec baffle_shape{ + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002, -BW), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002, 0.026 + Gate_height), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002 + 0.006, 0.026 + Gate_height), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002 + 0.006, -BW), + Vecd(Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002, -BW) +}; +/* Case-dependent geometries.*/ +class WaterBlock : public MultiPolygonShape +{ +public: + explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(water_block_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(baffle_shape, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(baffle_base_shap, ShapeBooleanOps::sub); + } +}; +/* Case-dependent geometries.*/ +class WallBoundary : public MultiPolygonShape +{ +public: + explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(baffle_base_shap, ShapeBooleanOps::add); + } +}; +/** Particle generator and constraint boundary for shell baffle. */ +class ShellBaffleParticleGenerator : public SurfaceParticleGenerator +{ +public: + explicit ShellBaffleParticleGenerator(SPHBody &sph_body) : SurfaceParticleGenerator(sph_body){}; + virtual void initializeGeometricVariables() override + { + for (int i = 0; i < particle_number_mid_surface; i++) + { + Real x = Gate_x - 0.0262 + 0.0075 + 0.008 + 0.006 + 0.002 + 0.006 - 0.4 * particle_spacing_ref; //0.5, 0.75, 1.5 + Real y = -BW + Real(i) * particle_spacing_ref; + initializePositionAndVolumetricMeasure(Vecd(x, y), particle_spacing_ref); + Vec2d normal_direction = Vec2d(1.0, 0); + initializeSurfaceProperties(normal_direction, thickness); + } + } +}; +/** Define the boundary geometry. */ +class BoundaryGeometry : public BodyPartByParticle +{ +public: + BoundaryGeometry(SPHBody &body, const std::string &body_part_name) + : BodyPartByParticle(body, body_part_name) + { + TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometry::tagManually, this, _1); + tagParticles(tagging_particle_method); + }; + virtual ~BoundaryGeometry(){}; + +private: + void tagManually(size_t index_i) + { + if (base_particles_.pos_[index_i][1] < 0.026) + { + body_part_particles_.push_back(index_i); + } + }; +}; +/** + * @brief define external force. + */ +class VariableGravity : public Gravity +{ +public: + VariableGravity() : Gravity(Vecd(0.0, -gravity_g)) {}; + void UpdateAcceleration() + { + + if (GlobalStaticVariables::physical_time_ <= 1.0) + { + global_acceleration_[0] = 0.0; + } + + if (GlobalStaticVariables::physical_time_ > 1.0) + { + global_acceleration_[0] + = k_0 * cos(w_0 * (GlobalStaticVariables::physical_time_ - 1.0)); + } + + } +}; +/** + * @brief Setup for gate and fluid boservers. + */ +Real h = 1.37 * particle_spacing_ref; +Real probe_x1 = DL - 0.032; +MultiPolygon CreateWaveProbeShape1() +{ + std::vector pnts; + pnts.push_back(Vecd(probe_x1 - h, 0.0)); + pnts.push_back(Vecd(probe_x1 - h, DH)); + pnts.push_back(Vecd(probe_x1 + h, DH)); + pnts.push_back(Vecd(probe_x1 + h, 0.0)); + pnts.push_back(Vecd(probe_x1 - h, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; +} + +Real probe_x2 = 0.032; +MultiPolygon CreateWaveProbeShape2() +{ + std::vector pnts; + pnts.push_back(Vecd(probe_x2 - h, 0.0)); + pnts.push_back(Vecd(probe_x2 - h, DH)); + pnts.push_back(Vecd(probe_x2 + h, DH)); + pnts.push_back(Vecd(probe_x2 + h, 0.0)); + pnts.push_back(Vecd(probe_x2 - h, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; +} + +Real probe_x3 = DL / 2.0; +MultiPolygon CreateWaveProbeShape3() +{ + std::vector pnts; + pnts.push_back(Vecd(probe_x3 - h, 0.0)); + pnts.push_back(Vecd(probe_x3 - h, DH)); + pnts.push_back(Vecd(probe_x3 + h, DH)); + pnts.push_back(Vecd(probe_x3 + h, 0.0)); + pnts.push_back(Vecd(probe_x3 - h, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; +} +/** + * @brief Define the observer body. + */ +StdVec baffle_disp_probe_location{Vecd(Gate_x, 0.195), Vecd(Gate_x, 0.145), Vecd(Gate_x, 0.095), Vecd(Gate_x, 0.05)}; +StdVec baffle_pressure_probe_location{Vecd(Gate_x - Gate_width / 2, 0.19), Vecd(Gate_x + Gate_width / 2, 0.19)}; +StdVec fluid_pressure_probe_location{Vecd(DL, 0.045), Vecd(DL, 0.095), Vecd(DL, 0.165), Vecd(DL, 0.22)}; \ No newline at end of file diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position.xml b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position.xml new file mode 100644 index 0000000000..1dbb9f836b --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position.xml @@ -0,0 +1,8732 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_0_result.xml b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_0_result.xml new file mode 100644 index 0000000000..ee29cf5e9c --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_0_result.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_5_result.xml b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_5_result.xml new file mode 100644 index 0000000000..4d96055559 --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_Run_5_result.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_dtwdistance.xml b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_dtwdistance.xml new file mode 100644 index 0000000000..3160c7e8d8 --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_runtimes.dat b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_runtimes.dat new file mode 100644 index 0000000000..de07de18db --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/BaffleDispObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +4 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/regression_test_tool.py b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..0bfc8bad6a --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/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_sloshing_baffle_fluid_shell +""" + +case_name = "test_2d_sloshing_baffle_fluid_shell" +body_name = "Observer" +parameter_name = "Position" + +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/2d_examples/test_2d_sloshing_baffle_fluid_shell/sloshing_baffle.cpp b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/sloshing_baffle.cpp new file mode 100644 index 0000000000..3110197a4c --- /dev/null +++ b/tests/2d_examples/test_2d_sloshing_baffle_fluid_shell/sloshing_baffle.cpp @@ -0,0 +1,228 @@ +/** +* @file sloshing-baffle.cpp +* @brief Fluid-shell interaction in sloshing flow. +* @details Here, the first fluid-shell interaction test is presented. +* @author Chi Zhang +*/ +#include "sphinxsys.h" +#include "case.h" +using namespace SPH; + +int main(int ac, char *av[]) +{ + /** + * @brief Build up -- a SPHSystem -- + */ + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + IOEnvironment io_environment(sph_system); +#ifdef BOOST_AVAILABLE + sph_system.handleCommandlineOptions(ac, av); +#endif + /** Set the starting time to zero. */ + GlobalStaticVariables::physical_time_ = 0.0; + /** The water block, body, material and particles container. */ + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); + water_block.generateParticles(); + /** The wall boundary, body and particles container. */ + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + wall_boundary.addBodyStateForRecording("NormalDirection"); + /** The baffle, body and particles container. */ + SolidBody shell_baffle(sph_system, makeShared("ShellBaffle")); + shell_baffle.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); + shell_baffle.generateParticles(); + shell_baffle.addBodyStateForRecording("PseudoNormal"); + /** @brief Particle and body creation of baffle observer.*/ + ObserverBody baffle_disp_observer(sph_system, "BaffleDispObserver"); + baffle_disp_observer.generateParticles(baffle_disp_probe_location); + /** Pressure probe on Flap. */ + ObserverBody baffle_pressure_observer(sph_system, "BafflePressureObserver"); + baffle_pressure_observer.generateParticles(baffle_pressure_probe_location); + /** @brief Particle and body creation of fluid observer.*/ + ObserverBody fluid_observer(sph_system, "Fluidobserver"); + fluid_observer.generateParticles(fluid_pressure_probe_location); + /** topology */ + InnerRelation water_block_inner(water_block); + InnerRelation baffle_inner(shell_baffle); + ContactRelation water_shell_contact(water_block, {&shell_baffle}); + ContactRelation water_wall_contact(water_block, {&wall_boundary}); + + ComplexRelation water_wall_complex(water_block_inner, water_wall_contact); + ComplexRelation water_shell_complex(water_block_inner, water_shell_contact); + + ContactRelation baffle_water_contact(shell_baffle, { &water_block }); + ContactRelation observer_contact_with_water(fluid_observer, { &water_block }); + ContactRelation observer_contact_with_baffle(baffle_disp_observer, { &shell_baffle}); + /** Density and wall norm. */ + SharedPtr gravity_ptr = makeShared(); + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + /** Algorithm for Fluid dynamics. */ + SimpleDynamics fluid_step_initialization(water_block, gravity_ptr); + ReduceDynamics fluid_advection_time_step(water_block, U_max); + ReduceDynamics fluid_acoustic_time_step(water_block); + InteractionWithUpdate update_fluid_density_by_summation(water_shell_contact, water_wall_complex); + Dynamics1Level fluid_pressure_relaxation(water_wall_contact, water_shell_complex); + Dynamics1Level fluid_density_relaxation(water_wall_contact, water_shell_complex); + /** Algorithms for solid. */ + ReduceDynamics shell_time_step_size(shell_baffle); + InteractionDynamics shell_corrected_configuration(baffle_inner); + Dynamics1Level shell_stress_relaxation_first(baffle_inner); + Dynamics1Level shell_stress_relaxation_second(baffle_inner); + /** FSI */ + InteractionDynamics viscous_force_on_shell(baffle_water_contact); + InteractionDynamics fluid_force_on_shell_update(baffle_water_contact, viscous_force_on_shell); + solid_dynamics::AverageVelocityAndAcceleration average_velocity_and_acceleration(shell_baffle); + /** constraint and damping */ + BoundaryGeometry shell_boundary_geometry(shell_baffle, "BoundaryGeometry"); + SimpleDynamics baffle_constrain(shell_boundary_geometry); + DampingWithRandomChoice>> + baffle_position_damping(0.2, baffle_inner, "Velocity", physical_viscosity); + DampingWithRandomChoice>> + baffle_rotation_damping(0.2, baffle_inner, "AngularVelocity", physical_viscosity); + /** + * @brief Output. + */ + BodyStatesRecordingToPlt write_real_body_states_to_plt(io_environment, sph_system.real_bodies_); + BodyStatesRecordingToVtp write_real_body_states_to_vtp(io_environment, sph_system.real_bodies_); + /** Output the observed displacement of baffle. */ + RegressionTestDynamicTimeWarping> + write_baffle_displacement("Position", io_environment, observer_contact_with_baffle); + /** Output the observed pressure of fluid. */ + RegressionTestDynamicTimeWarping> + write_fluid_pressure_wall("Pressure", io_environment, observer_contact_with_water); + /** WaveProbes. */ + /** #1.*/ + BodyRegionByCell wave_probe_buffer_no_1(water_block, makeShared(CreateWaveProbeShape1(), "WaveProbe_01")); + ReducedQuantityRecording> + wave_probe_1(io_environment, wave_probe_buffer_no_1); + /** #2.*/ + BodyRegionByCell wave_probe_buffer_no_2(water_block, makeShared(CreateWaveProbeShape2(), "WaveProbe_02")); + ReducedQuantityRecording> + wave_probe_2(io_environment, wave_probe_buffer_no_2); + /** #3.*/ + BodyRegionByCell wave_probe_buffer_no_3(water_block, makeShared(CreateWaveProbeShape3(), "WaveProbe_03")); + ReducedQuantityRecording> + wave_probe_3(io_environment, wave_probe_buffer_no_3); + /** + * @brief The time stepping starts here. + */ + /** Prepare quantities will be used once only and initial condition.*/ + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + wall_boundary_normal_direction.parallel_exec(); + shell_corrected_configuration.parallel_exec(); + /** Initial output. */ + //write_real_body_states_to_plt.writeToFile(0); + write_real_body_states_to_vtp.writeToFile(0); + wave_probe_1.writeToFile(0); + wave_probe_2.writeToFile(0); + wave_probe_3.writeToFile(0); + write_baffle_displacement.writeToFile(0); + write_fluid_pressure_wall.writeToFile(0); + /** Time parameters. */ + int number_of_iterations = 0; + int screen_output_interval = 100; + int observation_sample_interval = screen_output_interval * 2; + int restart_output_interval = screen_output_interval * 10; + Real End_Time = 5.0; /**< End time. */ + Real D_Time = End_Time /100; /**< time stamps for output. */ + Real Dt = 0.0; /**< Default advection time step sizes. */ + Real dt = 0.0; /**< Default acoustic time step sizes. */ + Real dt_s = 0.0; /**< Default acoustic time step sizes for solid. */ + tick_count t1 = tick_count::now(); + tick_count::interval_t interval; + /** + * @brief 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) + { + /** Acceleration due to viscous force and gravity. */ + gravity_ptr->UpdateAcceleration(); + fluid_step_initialization.parallel_exec(); + update_fluid_density_by_summation.parallel_exec(); + + viscous_force_on_shell.parallel_exec(); + + Dt = 0.5 * fluid_advection_time_step.parallel_exec(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + /** Fluid relaxation and force computaton. */ + fluid_pressure_relaxation.parallel_exec(dt); + fluid_force_on_shell_update.parallel_exec(); + fluid_density_relaxation.parallel_exec(dt); + + /** Solid dynamics time stepping. */ + Real dt_s_sum = 0.0; + average_velocity_and_acceleration.initialize_displacement_.parallel_exec(); + while (dt_s_sum < dt) + { + dt_s = shell_time_step_size.parallel_exec(); + if (dt - dt_s_sum < dt_s) dt_s = dt - dt_s_sum; + shell_stress_relaxation_first.parallel_exec(dt_s); + baffle_constrain.parallel_exec(); + shell_stress_relaxation_second.parallel_exec(dt_s); + dt_s_sum += dt_s; + } + average_velocity_and_acceleration.update_averages_.parallel_exec(dt); + + dt = fluid_acoustic_time_step.parallel_exec(); + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + } + + if (number_of_iterations % screen_output_interval == 0) + { + cout << fixed << setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " dt = " << dt << " dt_s = " << dt_s << "\n"; + } + number_of_iterations++; + + /** Update cell linked list and configuration. */ + water_block.updateCellLinkedList(); + shell_baffle.updateCellLinkedList(); + + water_block_inner.updateConfiguration(); + water_wall_contact.updateConfiguration(); + water_shell_contact.updateConfiguration(); + baffle_water_contact.updateConfiguration(); + observer_contact_with_water.updateConfiguration(); + + wave_probe_1.writeToFile(number_of_iterations); + wave_probe_2.writeToFile(number_of_iterations); + wave_probe_3.writeToFile(number_of_iterations); + write_baffle_displacement.writeToFile(number_of_iterations); + write_fluid_pressure_wall.writeToFile(number_of_iterations); + + } + write_real_body_states_to_vtp.writeToFile(); + + tick_count t2 = tick_count::now(); + + tick_count t3 = tick_count::now(); + interval += t3 - t2; + } + tick_count t4 = tick_count::now(); + tick_count::interval_t tt; + tt = t4 - t1 - interval; + cout << "Total wall time for computation: " << tt.seconds() << " seconds." << endl; + + if (sph_system.generate_regression_data_) + { + write_baffle_displacement.generateDataBase(1.0e-2); + } + else + { + write_baffle_displacement.newResultTest(); + } + + return 0; +} diff --git a/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp b/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp index b9a584bc4a..dd9b835d73 100644 --- a/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp +++ b/tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp @@ -76,7 +76,7 @@ int main() InteractionDynamics air_viscous_acceleration(air_water_complex); InteractionDynamics water_viscous_acceleration(water_air_complex); /** Surface tension and wetting effects. */ - InteractionWithUpdate surface_detection(water_air_complex.getInnerRelation(), water_wall_contact); + InteractionWithUpdate surface_detection(water_wall_contact, water_air_complex.getInnerRelation()); InteractionDynamics color_gradient(water_air_complex.getInnerRelation(), water_wall_contact); InteractionDynamics color_gradient_interpolation(water_air_complex.getInnerRelation()); InteractionDynamics surface_tension_acceleration(water_air_complex.getInnerRelation(), tension_force); diff --git a/tests/3d_examples/test_3d_poiseuille_flow/CMakeLists.txt b/tests/3d_examples/test_3d_poiseuille_flow/CMakeLists.txt new file mode 100644 index 0000000000..ccfa6747a9 --- /dev/null +++ b/tests/3d_examples/test_3d_poiseuille_flow/CMakeLists.txt @@ -0,0 +1,20 @@ +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}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_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}) + diff --git a/tests/3d_examples/test_3d_poiseuille_flow/poiseuille_flow.cpp b/tests/3d_examples/test_3d_poiseuille_flow/poiseuille_flow.cpp new file mode 100644 index 0000000000..159ff78484 --- /dev/null +++ b/tests/3d_examples/test_3d_poiseuille_flow/poiseuille_flow.cpp @@ -0,0 +1,348 @@ +/** + * @file poiseuille_flow.cpp + * @brief 3D poiseuille flow example + * @details This is the one of the basic test cases for validating viscous flow, with using emitter as inlet and disposer as outlet. + // TODO: include shell in next patch. + * @author + */ +/** + * @brief SPHinXsys Library. + */ +#include + +#include "base_data_type.h" +#include "sphinxsys.h" + +using namespace SPH; + + +class ObersverAxial : public ObserverParticleGenerator +{ +public: + ObersverAxial(SPHBody &sph_body, double full_length, Vec3d translation = Vec3d(0.0,0.0,0.0)) : ObserverParticleGenerator(sph_body) + { + int ny = 51; + for (int i = 0; i < ny; i++) + { + double y = full_length / (ny - 1) * i; + Vec3d point_coordinate(0.0, y, 0.0); + positions_.emplace_back(point_coordinate + translation); + } + } +}; + +class ObserverRadial : public ObserverParticleGenerator +{ +public: + ObserverRadial(SPHBody &sph_body, double full_length, double diameter, int number_of_particles, Vec3d translation = Vec3d(0.0,0.0,0.0)) : ObserverParticleGenerator(sph_body) + { + + int n = number_of_particles + 1; + double y = full_length / 2.0; + for (int i = 0; i < n - 1; i++) // we leave out the point clsoe to the boundary as the interpolation there is incorrect + { + double z = diameter / 2.0 * i / double(n); + positions_.emplace_back(Vec3d(0.0, y, z) + translation); + positions_.emplace_back(Vec3d(0.0, y, -z) + translation); + } + } +}; + +/** + * @brief Namespace cite here. + */ +using namespace SPH; +/** + * @brief Basic geometry parameters and numerical setup. + */ + +const Real scale = 0.001; +const Real diameter = 6.35 * scale; +const Real fluid_radius = 0.5 * diameter; +const Real full_length = 15 * fluid_radius; +const int number_of_particles = 10; +const Real resolution_ref = diameter / number_of_particles; +const Real inflow_length = resolution_ref * 20.0; // Inflow region +const Real wall_thickness = resolution_ref * 4.0; +const int simtk_resolution = 20; +const Vec3d translation_fluid(0.,full_length*0.5,0.); +/** + * @brief Geometry parameters for boundary condition. + */ +const Vec3d emitter_halfsize = Vec3d(fluid_radius, resolution_ref*2 ,fluid_radius); +const Vec3d emitter_translation = Vec3d(0., resolution_ref*2 ,0.); +const Vec3d emitter_buffer_halfsize = Vec3d(fluid_radius, resolution_ref*5 ,fluid_radius); +const Vec3d emitter_buffer_translation = Vec3d(0., resolution_ref*5 ,0.); +const Vec3d disposer_halfsize = Vec3d(fluid_radius*2, resolution_ref*2 ,fluid_radius*2); +const Vec3d disposer_translation = Vec3d(0., full_length, 0.) - Vec3d(0., disposer_halfsize[1], 0.); + +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec3d(-diameter , 0, -diameter) - Vec3d(wall_thickness, wall_thickness, wall_thickness), Vec3d(diameter , full_length, diameter) + Vec3d(wall_thickness, wall_thickness, wall_thickness)); +/** + * @brief Material properties of the fluid. + */ +const Real rho0_f = 1000.0; /**< Reference density of fluid. */ +const Real mu_f = 6.5e-3; /**< Viscosity. */ +const Real U_f = 30.0e-6 / 60. / (M_PI / 4.0 * diameter * diameter); /**< Characteristic velocity. Average velocity */ +const Real U_max = 2.0 * U_f; //parabolic inflow, Thus U_max = 2*U_f +const Real c_f = 10.0 * U_max; /**< Reference sound speed. */ +/** + * @brief Define water shape + */ +class WaterBlock : public ComplexShape +{ +public: + explicit WaterBlock(const std::string &shape_name) : ComplexShape(shape_name) + { + add(SimTK::UnitVec3(0., 1., 0.), fluid_radius, full_length*0.5, simtk_resolution, translation_fluid); + } +}; +/** + * @brief Define wall shape + */ +class WallBoundary : public ComplexShape +{ +public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + add(SimTK::UnitVec3(0.,1.,0.), fluid_radius + wall_thickness, full_length*0.5 + wall_thickness, simtk_resolution, translation_fluid); + subtract(SimTK::UnitVec3(0.,1.,0.), fluid_radius , full_length*0.5 + wall_thickness, simtk_resolution, translation_fluid); + } +}; +/** + * @brief Inflow velocity + */ +struct InflowVelocity +{ + Real u_ref_, t_ref_; + AlignedBoxShape &aligned_box_; + Vecd halfsize_; + + template + InflowVelocity(BoundaryConditionType &boundary_condition) + : u_ref_(U_f), t_ref_(2.0), + aligned_box_(boundary_condition.getAlignedBox()), + halfsize_(aligned_box_.HalfSize()) {} + + Vecd operator()(Vecd &position, Vecd &velocity) + { + Vecd target_velocity = Vec3d(0,0,0); + if (aligned_box_.checkInBounds(0, position)) + { + target_velocity[1] = 2.0 * U_f * (1.0 - (position[0] * position[0] + position[2] * position[2]) / fluid_radius/ fluid_radius); + } + return target_velocity; + } +}; + +/** + * @brief Main program starts here. + */ +int main() +{ + /** + * @brief Build up -- a SPHSystem -- + */ + SPHSystem system(system_domain_bounds, resolution_ref); + /** Set the starting time. */ + GlobalStaticVariables::physical_time_ = 0.0; + /** + * @brief Material property, particles and body creation of fluid. + */ + FluidBody water_block(system, makeShared("WaterBody")); + water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); + water_block.generateParticles(); + /** + * @brief Particle and body creation of wall boundary. + */ + SolidBody wall_boundary(system, makeShared("Wall")); + wall_boundary.defineAdaptation(1.15, 1.0); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + /** topology */ + ComplexRelation water_block_complex(water_block, {&wall_boundary}); + /** + * @brief Define all numerical methods which are used in this case. + */ + /** + * @brief Methods used for time stepping. + */ + /** Define external force. */ + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + /** Initialize particle acceleration. */ + SimpleDynamics initialize_a_fluid_step(water_block, makeShared(Vec3d(0.0, 0.0, 0.0))); + /** + * @brief Algorithms of fluid dynamics. + */ + /** time-space method to detect surface particles. (Important for DensitySummationFreeSurfaceComplex work correctly.)*/ + InteractionWithUpdate inlet_outlet_surface_particle_indicator(water_block_complex); + /** Evaluation of density by summation approach. */ + InteractionWithUpdate update_density_by_summation(water_block_complex); + /** Time step size without considering sound wave speed. */ + ReduceDynamics get_fluid_advection_time_step_size(water_block, U_max); + /** Time step size with considering sound wave speed. */ + ReduceDynamics get_fluid_time_step_size(water_block); + /** Pressure relaxation algorithm without Riemann solver for viscous flows. */ + Dynamics1Level pressure_relaxation(water_block_complex); + /** Pressure relaxation algorithm by using position verlet time stepping. */ + Dynamics1Level density_relaxation(water_block_complex); + /** Computing viscous acceleration. */ + InteractionDynamics viscous_acceleration(water_block_complex); + /** Impose transport velocity. */ + InteractionDynamics transport_velocity_correction(water_block_complex); + /** + * @brief Boundary conditions. Inflow & Outflow in Y-direction + */ + BodyAlignedBoxByParticle emitter( + water_block, makeShared(Transform3d(Vec3d(emitter_translation)), emitter_halfsize)); + SimpleDynamics emitter_inflow_injection(emitter, 10, 1); + /** Emitter buffer inflow condition. */ + BodyAlignedBoxByCell emitter_buffer( + water_block, makeShared(Transform3d(Vec3d(emitter_buffer_translation)), emitter_buffer_halfsize)); + SimpleDynamics, BodyAlignedBoxByCell> emitter_buffer_inflow_condition(emitter_buffer); + BodyAlignedBoxByCell disposer( + water_block, makeShared(Transform3d(Vec3d(disposer_translation)), disposer_halfsize)); + SimpleDynamics disposer_outflow_deletion(disposer, 1); + + /** + * @brief Output. + */ + IOEnvironment io_environment(system); + water_block.addBodyStateForRecording("SurfaceIndicator"); + /** Output the body states. */ + BodyStatesRecordingToVtp body_states_recording(io_environment, system.real_bodies_); + /** + * @brief OBSERVER. + */ + ObserverBody observer_axial(system, "fluid_observer_axial"); + observer_axial.generateParticles(full_length); + ObserverBody observer_radial(system, "fluid_observer_radial"); + observer_radial.generateParticles(full_length, diameter, number_of_particles); + ContactRelation observer_contact_axial(observer_axial, {&water_block}); + ContactRelation observer_contact_radial(observer_radial, {&water_block}); + ObservedQuantityRecording write_fluid_velocity_axial("Velocity", io_environment, observer_contact_axial); + ObservedQuantityRecording write_fluid_velocity_radial("Velocity", io_environment, observer_contact_radial); + /** + * @brief Setup geometry and initial conditions. + */ + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + wall_boundary_normal_direction.parallel_exec(); + /** Output the start states of bodies. */ + body_states_recording.writeToFile(0); + /** + * @brief Basic parameters. + */ + size_t number_of_iterations = system.RestartStep(); + int screen_output_interval = 100; + Real end_time = 10.0; /**< End time. */ + Real Output_Time = 0.1; /**< Time stamps for output of body states. */ + Real dt = 0.0; /**< Default acoustic time step sizes. */ + /** statistics for computing CPU time. */ + tick_count t1 = tick_count::now(); + tick_count::interval_t interval; + tick_count::interval_t interval_computing_time_step; + tick_count::interval_t interval_computing_pressure_relaxation; + tick_count::interval_t interval_updating_configuration; + tick_count time_instance; + /** + * @brief 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 < Output_Time) + { + /** Acceleration due to viscous force and gravity. */ + time_instance = tick_count::now(); + initialize_a_fluid_step.parallel_exec(); + Real Dt = get_fluid_advection_time_step_size.parallel_exec(); + inlet_outlet_surface_particle_indicator.parallel_exec(); + update_density_by_summation.parallel_exec(); + viscous_acceleration.parallel_exec(); + transport_velocity_correction.parallel_exec(); + interval_computing_time_step += tick_count::now() - time_instance; + /** Dynamics including pressure relaxation. */ + time_instance = tick_count::now(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + dt = SMIN(get_fluid_time_step_size.parallel_exec(), Dt - relaxation_time); + pressure_relaxation.parallel_exec(dt); + emitter_buffer_inflow_condition.parallel_exec(); + density_relaxation.parallel_exec(dt); + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + } + interval_computing_pressure_relaxation += tick_count::now() - time_instance; + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " dt = " << dt << "\n"; + body_states_recording.writeToFile(); + + } + number_of_iterations++; + /** Update cell linked list and configuration. */ + time_instance = tick_count::now(); + + /** Water block configuration and periodic condition. */ + emitter_inflow_injection.parallel_exec(); + disposer_outflow_deletion.parallel_exec(); + + + water_block.updateCellLinkedListWithParticleSort(100); + water_block_complex.updateConfiguration(); + interval_updating_configuration += tick_count::now() - time_instance; + } + tick_count t2 = tick_count::now(); + tick_count t3 = tick_count::now(); + interval += t3 - t2; + + /** Update observer and write output of observer. */ + observer_contact_axial.updateConfiguration(); + observer_contact_radial.updateConfiguration(); + write_fluid_velocity_axial.writeToFile(number_of_iterations); + write_fluid_velocity_radial.writeToFile(number_of_iterations); + } + tick_count t4 = tick_count::now(); + + tick_count::interval_t tt; + tt = t4 - t1 - interval; + std::cout << "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_pressure_relaxation = " + << interval_computing_pressure_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + /** + * @brief Gtest strat from here. + */ + /* Define analytical solution of the inflow velocity.*/ + std::function inflow_velocity = [&](Vec3d pos){ + return Vec3d( + 0.0, + 2.0 * U_f * (1.0 - (pos[0] * pos[0] + pos[2] * pos[2]) / (diameter*0.5) / (diameter*0.5)), + 0.0 + ); + }; + /* Compare all simulation to the anayltical solution. */ + //Axial direction. + for (size_t i = 0; i < observer_axial.getBaseParticles().pos_.size(); i++) + { + EXPECT_NEAR(inflow_velocity(observer_axial.getBaseParticles().pos_[i])[1], observer_axial.getBaseParticles().vel_[i][1], U_max * 10e-2); // it's below 5% but 10% for CI + } + // Radial direction + for (size_t i = 0; i < observer_radial.getBaseParticles().pos_.size(); i++) + { + EXPECT_NEAR(inflow_velocity(observer_radial.getBaseParticles().pos_[i])[1], observer_radial.getBaseParticles().vel_[i][1], U_max * 10e-2); // it's below 5% but 10% for CI + } + + + return 0; +} diff --git a/tests/3d_examples/test_3d_poiseuille_flow_shell/CMakeLists.txt b/tests/3d_examples/test_3d_poiseuille_flow_shell/CMakeLists.txt new file mode 100644 index 0000000000..ccfa6747a9 --- /dev/null +++ b/tests/3d_examples/test_3d_poiseuille_flow_shell/CMakeLists.txt @@ -0,0 +1,20 @@ +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}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_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}) + diff --git a/tests/3d_examples/test_3d_poiseuille_flow_shell/poiseuille_flow_shell.cpp b/tests/3d_examples/test_3d_poiseuille_flow_shell/poiseuille_flow_shell.cpp new file mode 100644 index 0000000000..2fe2e16d1e --- /dev/null +++ b/tests/3d_examples/test_3d_poiseuille_flow_shell/poiseuille_flow_shell.cpp @@ -0,0 +1,423 @@ +/** + * @file poiseuille_flow.cpp + * @brief 3D poiseuille flow example + * @details This is the one of the basic test cases for validating viscous flow, + with using emitter as inlet and disposer as outlet. + // TODO: include shell in next patch. + * @author + */ +/** + * @brief SPHinXsys Library. + */ +#include "base_data_type.h" +#include "sphinxsys.h" +#include +using namespace SPH; + +class ObersverAxial : public ObserverParticleGenerator { +public: + ObersverAxial(SPHBody &sph_body, double full_length, + Vec3d translation = Vec3d(0.0, 0.0, 0.0)) + : ObserverParticleGenerator(sph_body) { + int ny = 51; + for (int i = 0; i < ny; i++) { + double y = full_length / (ny - 1) * i; + Vec3d point_coordinate(0.0, y, 0.0); + positions_.emplace_back(point_coordinate + translation); + } + } +}; + +class ObserverRadial : public ObserverParticleGenerator { +public: + ObserverRadial(SPHBody &sph_body, double full_length, double diameter, + int number_of_particles, + Vec3d translation = Vec3d(0.0, 0.0, 0.0)) + : ObserverParticleGenerator(sph_body) { + + int n = number_of_particles + 1; + double y = full_length / 2.0; + for (int i = 0; i < n - 1; + i++) // we leave out the point clsoe to the boundary as the + // interpolation there is incorrect + { + double z = diameter / 2.0 * i / double(n); + positions_.emplace_back(Vec3d(0.0, y, z) + translation); + positions_.emplace_back(Vec3d(0.0, y, -z) + translation); + } + } +}; + +/** + * @brief Namespace cite here. + */ +using namespace SPH; +/** + * @brief Basic geometry parameters and numerical setup. + */ + +const Real scale = 0.001; +const Real diameter = 6.35 * scale * 1; +const Real fluid_radius = 0.5 * diameter; +const Real full_length = 15 * fluid_radius; +const int number_of_particles = 10; +const Real resolution_ref = diameter / number_of_particles; +const Real inflow_length = resolution_ref * 20.0; // Inflow region +const Real wall_thickness = resolution_ref * 4.0; +const Real shell_thickness = resolution_ref * 2.0; +const int simtk_resolution = 20; +const Vec3d translation_fluid(0., full_length * 0.5, 0.); +/** + * @brief Geometry parameters for shell. + */ +const Real radius_mid_surface = fluid_radius + shell_thickness * 0.5; +const int particle_number_mid_surface = + int(2.0 * radius_mid_surface * Pi / resolution_ref); +const int particle_number_height = + 2 * int((full_length * 0.5 + wall_thickness) / resolution_ref); +const int BWD = + 1; /** Width of the boundary layer measured by number of particles. */ + +/** + * @brief Geometry parameters for boundary condition. + */ +const Vec3d emitter_halfsize = + Vec3d(fluid_radius, resolution_ref * 2, fluid_radius); +const Vec3d emitter_translation = Vec3d(0., resolution_ref * 2, 0.); +const Vec3d emitter_buffer_halfsize = + Vec3d(fluid_radius, inflow_length * 0.5, fluid_radius); +const Vec3d emitter_buffer_translation = Vec3d(0., inflow_length * 0.5, 0.); +const Vec3d disposer_halfsize = + Vec3d(fluid_radius * 2, resolution_ref * 2, fluid_radius * 2); +const Vec3d disposer_translation = + Vec3d(0., full_length, 0.) - Vec3d(0., disposer_halfsize[1], 0.); + +/** Domain bounds of the system. */ +const BoundingBox system_domain_bounds(Vec3d(-diameter, 0, -diameter) - + Vec3d(wall_thickness, wall_thickness, + wall_thickness), + Vec3d(diameter, full_length, diameter) + + Vec3d(wall_thickness, wall_thickness, + wall_thickness)); +/** + * @brief Material properties of the fluid. + */ +const Real rho0_f = 1000.0; /**< Reference density of fluid. */ +const Real mu_f = 6.5e-3; /**< Viscosity. */ +const Real Re = 10; +// const Real U_f = 30.0e-6 / 60. / (M_PI / 4.0 * diameter * diameter); /**< +// Characteristic velocity. Average velocity */ +const Real U_f = Re / diameter * mu_f * 0.5; +const Real U_max = 2.0 * U_f; // parabolic inflow, Thus U_max = 2*U_f +const Real c_f = 10.0 * U_max; /**< Reference sound speed. */ +/** + * @brief Define water shape + */ +class WaterBlock : public ComplexShape { +public: + explicit WaterBlock(const std::string &shape_name) + : ComplexShape(shape_name) { + add(SimTK::UnitVec3(0., 1., 0.), fluid_radius, + full_length * 0.5, simtk_resolution, + translation_fluid); + } +}; +/** + * @brief Define wall shape + */ +class ShellBoundary : public SurfaceParticleGenerator { +public: + explicit ShellBoundary(SPHBody &sph_body) + : SurfaceParticleGenerator(sph_body){}; + void initializeGeometricVariables() override { + for (int i = 0; i < particle_number_mid_surface + 2 * BWD; i++) { + for (int j = 0; j < particle_number_height; j++) { + Real x = radius_mid_surface * + cos(Pi + (i - BWD + 0.5) * 2 * Pi / + (Real)particle_number_mid_surface); + Real y = (j - particle_number_height / 2) * resolution_ref + + resolution_ref * 0.5 + full_length * 0.5; + Real z = radius_mid_surface * + sin(Pi + (i - BWD + 0.5) * 2 * Pi / + (Real)particle_number_mid_surface); + initializePositionAndVolumetricMeasure(Vec3d(x, y, z), + resolution_ref * resolution_ref); + Vec3d n_0 = Vec3d(x / radius_mid_surface, 0.0, z / radius_mid_surface); + initializeSurfaceProperties(n_0, shell_thickness); + } + } + } +}; +/** + * @brief Inflow velocity + */ +struct InflowVelocity { + Real u_ref_, t_ref_; + AlignedBoxShape &aligned_box_; + Vec3d halfsize_; + + template + InflowVelocity(BoundaryConditionType &boundary_condition) + : u_ref_(U_f), t_ref_(2.0), + aligned_box_(boundary_condition.getAlignedBox()), + halfsize_(aligned_box_.HalfSize()) {} + + Vec3d operator()(Vec3d &position, Vec3d &velocity) { + Vec3d target_velocity = Vec3d(0, 0, 0); + if (aligned_box_.checkInBounds(0, position)) { + target_velocity[1] = + 2.0 * U_f * + (1.0 - (position[0] * position[0] + position[2] * position[2]) / + fluid_radius / fluid_radius); + } + return target_velocity; + } +}; + +/** + * @brief Main program starts here. + */ +int main() { + /** + * @brief Build up -- a SPHSystem -- + */ + SPHSystem system(system_domain_bounds, resolution_ref); + /** Set the starting time. */ + GlobalStaticVariables::physical_time_ = 0.0; + /** + * @brief Material property, particles and body creation of fluid. + */ + FluidBody water_block(system, makeShared("WaterBody")); + water_block + .defineParticlesAndMaterial( + rho0_f, c_f, mu_f); + water_block.generateParticles(); + /** + * @brief Particle and body creation of wall boundary. + */ + SolidBody shell_boundary(system, makeShared("Shell")); + shell_boundary.defineParticlesAndMaterial( + rho0_f, 1e3, 0.45); + shell_boundary.generateParticles(); + /** topology */ + InnerRelation water_inner(water_block); + ComplexRelation water_block_complex(water_block, {&shell_boundary}); + /** + * @brief Define all numerical methods which are used in this case. + */ + /** + * @brief Methods used for time stepping. + */ + /** Define external force. */ + SimpleDynamics shell_boundary_normal_direction( + shell_boundary); + /** Initialize particle acceleration. */ + SimpleDynamics initialize_a_fluid_step( + water_block, makeShared(Vec3d(0.0, 0.0, 0.0))); + /** + * @brief Algorithms of fluid dynamics. + */ + /** time-space method to detect surface particles. (Important for + * DensitySummationFreeSurfaceComplex work correctly.)*/ + InteractionWithUpdate< + fluid_dynamics::SpatialTemporalFreeSurfaceIdentificationComplex> + inlet_outlet_surface_particle_indicator(water_block_complex); + /** Evaluation of density by summation approach. */ + InteractionWithUpdate + update_density_by_summation(water_block_complex); + /** Time step size without considering sound wave speed. */ + ReduceDynamics + get_fluid_advection_time_step_size(water_block, U_max); + /** Time step size with considering sound wave speed. */ + ReduceDynamics get_fluid_time_step_size( + water_block); + /** Pressure relaxation algorithm without Riemann solver for viscous flows. */ + Dynamics1Level + pressure_relaxation(water_block_complex); + /** Pressure relaxation algorithm by using position verlet time stepping. */ + Dynamics1Level + density_relaxation(water_block_complex); + /** Computing viscous acceleration. */ + InteractionDynamics + viscous_acceleration(water_block_complex); + /** Impose transport velocity. */ + InteractionDynamics + transport_velocity_correction(water_block_complex); + /** + * @brief Boundary conditions. Inflow & Outflow in Y-direction + */ + BodyAlignedBoxByParticle emitter( + water_block, + makeShared(Transform3d(Vec3d(emitter_translation)), + emitter_halfsize)); + SimpleDynamics + emitter_inflow_injection(emitter, 10, 1); + /** Emitter buffer inflow condition. */ + BodyAlignedBoxByCell emitter_buffer( + water_block, makeShared( + Transform3d(Vec3d(emitter_buffer_translation)), + emitter_buffer_halfsize)); + SimpleDynamics, + BodyAlignedBoxByCell> + emitter_buffer_inflow_condition(emitter_buffer); + BodyAlignedBoxByCell disposer( + water_block, + makeShared(Transform3d(Vec3d(disposer_translation)), + disposer_halfsize)); + SimpleDynamics + disposer_outflow_deletion(disposer, 1); + + /** + * @brief Output. + */ + IOEnvironment io_environment(system); + water_block.addBodyStateForRecording("SurfaceIndicator"); + + shell_boundary.addBodyStateForRecording("MassiveMeasure"); + shell_boundary.addBodyStateForRecording("Density"); + shell_boundary.addBodyStateForRecording("VolumetricMeasure"); + shell_boundary.addBodyStateForRecording("Thickness"); + /** Output the body states. */ + BodyStatesRecordingToVtp body_states_recording(io_environment, + system.real_bodies_); + /** + * @brief OBSERVER. + */ + ObserverBody observer_axial(system, "fluid_observer_axial"); + observer_axial.generateParticles(full_length); + ObserverBody observer_radial(system, "fluid_observer_radial"); + observer_radial.generateParticles(full_length, diameter, + number_of_particles); + ContactRelation observer_contact_axial(observer_axial, {&water_block}); + ContactRelation observer_contact_radial(observer_radial, {&water_block}); + ObservedQuantityRecording write_fluid_velocity_axial( + "Velocity", io_environment, observer_contact_axial); + ObservedQuantityRecording write_fluid_velocity_radial( + "Velocity", io_environment, observer_contact_radial); + /** + * @brief Setup geometry and initial conditions. + */ + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + shell_boundary_normal_direction.parallel_exec(); + /** Output the start states of bodies. */ + body_states_recording.writeToFile(0); + /** + * @brief Basic parameters. + */ + size_t number_of_iterations = system.RestartStep(); + int screen_output_interval = 100; + Real end_time = 10.0; /**< End time. */ + Real Output_Time = 0.1; /**< Time stamps for output of body states. */ + Real dt = 0.0; /**< Default acoustic time step sizes. */ + /** statistics for computing CPU time. */ + tick_count t1 = tick_count::now(); + tick_count::interval_t interval; + tick_count::interval_t interval_computing_time_step; + tick_count::interval_t interval_computing_pressure_relaxation; + tick_count::interval_t interval_updating_configuration; + tick_count time_instance; + /** + * @brief 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 < Output_Time) { + /** Acceleration due to viscous force and gravity. */ + time_instance = tick_count::now(); + initialize_a_fluid_step.parallel_exec(); + Real Dt = get_fluid_advection_time_step_size.parallel_exec(); + inlet_outlet_surface_particle_indicator.parallel_exec(); + update_density_by_summation.parallel_exec(); + viscous_acceleration.parallel_exec(); + transport_velocity_correction.parallel_exec(); + interval_computing_time_step += tick_count::now() - time_instance; + /** Dynamics including pressure relaxation. */ + time_instance = tick_count::now(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) { + dt = SMIN(get_fluid_time_step_size.parallel_exec(), + Dt - relaxation_time); + pressure_relaxation.parallel_exec(dt); + density_relaxation.parallel_exec(dt); + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + emitter_buffer_inflow_condition.parallel_exec(); + } + interval_computing_pressure_relaxation += + tick_count::now() - time_instance; + if (number_of_iterations % screen_output_interval == 0) { + std::cout << std::fixed << std::setprecision(9) + << "N=" << number_of_iterations + << " Time = " << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " dt = " << dt << "\n"; + body_states_recording.writeToFile(); + } + number_of_iterations++; + /** Update cell linked list and configuration. */ + time_instance = tick_count::now(); + + /** Water block configuration and periodic condition. */ + emitter_inflow_injection.parallel_exec(); + disposer_outflow_deletion.parallel_exec(); + + water_block.updateCellLinkedListWithParticleSort(100); + water_block_complex.updateConfiguration(); + interval_updating_configuration += tick_count::now() - time_instance; + } + tick_count t2 = tick_count::now(); + tick_count t3 = tick_count::now(); + interval += t3 - t2; + + /** Update observer and write output of observer. */ + observer_contact_axial.updateConfiguration(); + observer_contact_radial.updateConfiguration(); + write_fluid_velocity_axial.writeToFile(number_of_iterations); + write_fluid_velocity_radial.writeToFile(number_of_iterations); + } + tick_count t4 = tick_count::now(); + + tick_count::interval_t tt; + tt = t4 - t1 - interval; + std::cout << "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_pressure_relaxation = " + << interval_computing_pressure_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) + << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + + /** + * @brief Gtest strat from here. + */ + /* Define analytical solution of the inflow velocity.*/ + std::function inflow_velocity = [&](Vec3d pos) { + return Vec3d(0.0, + 2.0 * U_f * + (1.0 - (pos[0] * pos[0] + pos[2] * pos[2]) / + (diameter * 0.5) / (diameter * 0.5)), + 0.0); + }; + /* Compare all simulation to the anayltical solution. */ + // Axial direction. + for (size_t i = 0; i < observer_axial.getBaseParticles().pos_.size(); i++) { + EXPECT_NEAR(inflow_velocity(observer_axial.getBaseParticles().pos_[i])[1], + observer_axial.getBaseParticles().vel_[i][1], + U_max * 10e-2); // it's below 5% but 10% for CI + } + // Radial direction + for (size_t i = 0; i < observer_radial.getBaseParticles().pos_.size(); i++) { + EXPECT_NEAR(inflow_velocity(observer_radial.getBaseParticles().pos_[i])[1], + observer_radial.getBaseParticles().vel_[i][1], + U_max * 10e-2); // it's below 5% but 10% for CI + } + + return 0; +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ab5f6d1fa2..99675d996e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,7 @@ option(SPHINXSYS_BUILD_2D_EXAMPLES "SPHINXSYS_BUILD_2D_EXAMPLES" ON) option(SPHINXSYS_BUILD_3D_EXAMPLES "SPHINXSYS_BUILD_3D_EXAMPLES" ON) option(SPHINXSYS_BUILD_UNIT_TESTS "SPHINXSYS_BUILD_UNIT_TESTS" ON) +option(SPHINXSYS_BUILD_USER_TESTS "SPHINXSYS_BUILD_USER_TESTS" ON) find_package(GTest CONFIG REQUIRED) include(GoogleTest) @@ -15,6 +16,10 @@ if(SPHINXSYS_BUILD_UNIT_TESTS) ADD_SUBDIRECTORY(unit_tests_src) endif() +if(SPHINXSYS_BUILD_USER_TESTS) + ADD_SUBDIRECTORY(user_examples) +endif() + add_subdirectory(modules) if(SPHINXSYS_3D AND SPHINXSYS_BUILD_3D_EXAMPLES) diff --git a/tests/user_examples/test_2d_square_fsi/2d_square_fsi.cpp b/tests/user_examples/test_2d_square_fsi/2d_square_fsi.cpp new file mode 100644 index 0000000000..201aa34f66 --- /dev/null +++ b/tests/user_examples/test_2d_square_fsi/2d_square_fsi.cpp @@ -0,0 +1,182 @@ +#include "sphinxsys.h" +#include "case.h" + +using namespace SPH; + +int main(int ac, char *av[]) +{ + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + /** handle command line arguments. */ + sph_system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(sph_system); + /** Set the starting time to zero. */ + GlobalStaticVariables::physical_time_ = 0.0; + /** The water block, body, material and particles container. */ + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); + water_block.generateParticles(); + water_block.addBodyStateForRecording("Pressure"); + water_block.addBodyStateForRecording("SurfaceIndicator"); + /** The wall boundary, body and particles container. */ + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + wall_boundary.addBodyStateForRecording("NormalDirection"); + /** The baffle, body and particles container. */ + SolidBody shell_baffle(sph_system, makeShared("ShellBaffle")); + shell_baffle.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); + shell_baffle.generateParticles(); + shell_baffle.addBodyStateForRecording("PseudoNormal"); + /** Particle and body creation of baffle observer.*/ + ObserverBody fluid_observer(sph_system, "FluidObserver"); + fluid_observer.generateParticles(observation_locations); + /** topology */ + InnerRelation water_inner(water_block); + InnerRelation baffle_inner(shell_baffle); + + ContactRelation water_wall_contact(water_block, {&wall_boundary}); + + ComplexRelation water_shell_complex(water_inner, {&shell_baffle}); + ContactRelation shell_water_contact(shell_baffle, {&water_block}); + + ContactRelation fluid_observer_contact(fluid_observer, {&water_block}); + /** Emitter setup. */ + BodyAlignedBoxByParticle emitter(water_block, makeShared(Transform2d(Vec2d(emitter_translation)), emitter_halfsize)); + SimpleDynamics emitter_inflow_injection(emitter, 10, 0); + /** Emitter buffer inflow condition. */ + BodyAlignedBoxByCell emitter_buffer(water_block, makeShared(Transform2d(Vec2d(emitter_buffer_translation)), emitter_buffer_halfsize)); + SimpleDynamics, BodyAlignedBoxByCell> emitter_buffer_inflow_condition(emitter_buffer); + + BodyAlignedBoxByCell disposer(water_block, makeShared(Transform2d(Vec2d(disposer_translation)), disposer_halfsize)); + SimpleDynamics disposer_outflow_deletion(disposer, 0); + /** Free-stream BC. */ + InteractionWithUpdate free_stream_surface_indicator(water_wall_contact, water_shell_complex); + InteractionWithUpdate update_fluid_density(water_wall_contact, water_shell_complex); + /** Algorithms for Fluid dynamics. */ + SimpleDynamics initialize_a_fluid_step(water_block, makeShared(Vec2d(0))); + ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); + ReduceDynamics get_fluid_time_step_size(water_block); + SimpleDynamics> velocity_boundary_condition_constraint(water_block); + Dynamics1Level pressure_relaxation(water_wall_contact, water_shell_complex); + pressure_relaxation.post_processes_.push_back(&velocity_boundary_condition_constraint); + Dynamics1Level density_relaxation(water_wall_contact, water_shell_complex); + InteractionDynamics viscous_acceleration(water_wall_contact, water_shell_complex); + InteractionDynamics transport_velocity_correction(water_wall_contact, water_shell_complex); + InteractionDynamics compute_vorticity(water_inner); + /** Algorithms for solid. */ + ReduceDynamics shell_time_step_size(shell_baffle); + InteractionDynamics shell_corrected_configuration(baffle_inner); + Dynamics1Level shell_stress_relaxation_first(baffle_inner, 5, true); + Dynamics1Level shell_stress_relaxation_second(baffle_inner); + /** Algorithms of FSI. */ + /** Compute the force exerted on solid body due to fluid pressure and viscosity. */ + InteractionDynamics viscous_force_on_shell(shell_water_contact); + InteractionDynamics fluid_force_on_shell_update(shell_water_contact, viscous_force_on_shell); + solid_dynamics::AverageVelocityAndAcceleration average_velocity_and_acceleration(shell_baffle); + /** constraint and damping */ + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + BoundaryGeometry shell_boundary_geometry(shell_baffle, "BoundaryGeometry"); + SimpleDynamics baffle_constrain(shell_boundary_geometry); + /** IO */ + BodyStatesRecordingToVtp write_real_body_states(io_environment, sph_system.real_bodies_); + RestartIO restart_io(io_environment, sph_system.real_bodies_); + ObservedQuantityRecording write_fluid_velocity("Velocity", io_environment, fluid_observer_contact); + /** 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.parallel_exec(); + shell_corrected_configuration.parallel_exec(); + write_real_body_states.writeToFile(); + /** Time parameters. */ + size_t number_of_iterations = 0; + int screen_output_interval = 100; + int restart_output_interval = screen_output_interval * 10; + Real end_time = 400.0; + Real output_interval = end_time / 200.0; + tick_count t1 = tick_count::now(); + tick_count::interval_t interval; + /** 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 < output_interval) + { + initialize_a_fluid_step.parallel_exec(); + Real Dt = get_fluid_advection_time_step_size.parallel_exec(); + free_stream_surface_indicator.parallel_exec(); + update_fluid_density.parallel_exec(); + viscous_acceleration.parallel_exec(); + transport_velocity_correction.parallel_exec(); + + /** FSI for viscous force. */ + viscous_force_on_shell.parallel_exec(); + + size_t inner_ite_dt = 0; + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + Real dt = SMIN(get_fluid_time_step_size.parallel_exec(), Dt - relaxation_time); + /** Fluid dynamics and force on solid. */ + pressure_relaxation.parallel_exec(dt); + fluid_force_on_shell_update.parallel_exec(); + density_relaxation.parallel_exec(dt); + + /** Solid dynamics time stepping. */ + Real dt_s_sum = 0.0; + average_velocity_and_acceleration.initialize_displacement_.parallel_exec(); + while (dt_s_sum < dt) + { + Real dt_s = shell_time_step_size.parallel_exec(); + if (dt - dt_s_sum < dt_s) dt_s = dt - dt_s_sum; + shell_stress_relaxation_first.parallel_exec(dt_s); + baffle_constrain.parallel_exec(); + shell_stress_relaxation_second.parallel_exec(dt_s); + dt_s_sum += dt_s; + } + average_velocity_and_acceleration.update_averages_.parallel_exec(dt); + + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + emitter_buffer_inflow_condition.parallel_exec(); + inner_ite_dt++; + } + + if (number_of_iterations % screen_output_interval == 0) + { + cout << fixed << setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " Dt / dt = " << inner_ite_dt << "\n"; + } + number_of_iterations++; + + /** Water block configuration and periodic condition. */ + emitter_inflow_injection.parallel_exec(); + disposer_outflow_deletion.parallel_exec(); + + water_block.updateCellLinkedListWithParticleSort(100); + shell_baffle.updateCellLinkedList(); + + water_shell_complex.updateConfiguration(); + water_wall_contact.updateConfiguration(); + shell_water_contact.updateConfiguration(); + } + + tick_count t2 = tick_count::now(); + /** write run-time observation into file */ + compute_vorticity.parallel_exec(); + write_real_body_states.writeToFile(); + fluid_observer_contact.updateConfiguration(); + write_fluid_velocity.writeToFile(number_of_iterations); + tick_count t3 = tick_count::now(); + interval += t3 - t2; + } + tick_count t4 = tick_count::now(); + + tick_count::interval_t tt; + tt = t4 - t1 - interval; + cout << "Total wall time for computation: " << tt.seconds() << " seconds." << endl; + + return 0; +} diff --git a/tests/user_examples/test_2d_square_fsi/CMakeLists.txt b/tests/user_examples/test_2d_square_fsi/CMakeLists.txt new file mode 100644 index 0000000000..e41c05c620 --- /dev/null +++ b/tests/user_examples/test_2d_square_fsi/CMakeLists.txt @@ -0,0 +1,18 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir + +set(CMAKE_VERBOSE_MAKEFILE on) + +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}) + +aux_source_directory(. DIR_SRCS) +add_executable(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS} ) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} sphinxsys_2d) \ No newline at end of file diff --git a/tests/user_examples/test_2d_square_fsi/case.h b/tests/user_examples/test_2d_square_fsi/case.h new file mode 100644 index 0000000000..dd75eabdde --- /dev/null +++ b/tests/user_examples/test_2d_square_fsi/case.h @@ -0,0 +1,166 @@ +/** + * @file case.h + * @brief This is the case study the flow induced vibration of a flexible thin structure + * attached to a square. + * @author Chi Zhang & Xiangyu Hu + */ + +#include "sphinxsys.h" +using namespace SPH; + +Real DL = 19.5; /**< Channel length. */ +Real DH = 8; /**< Channel height. */ +Vec2d insert_circle_center(5.0, 0.5 * DH); /**< Location of the cylinder center. */ +Real insert_circle_radius = 0.5; /**< Radius of the cylinder. */ +Real bh = 0.1 * insert_circle_radius; /**< Height of the beam. */ +Real bl = 8.0 * insert_circle_radius; /**< Length of the beam. */ +Real particle_spacing_ref = bh; /**< Global reference resolution. */ +Real DL_sponge = particle_spacing_ref * 20.0; /**< Sponge region to impose inflow condition. */ +Real BW = particle_spacing_ref * 4.0; /**< Boundary width, determined by specific layer of boundary particles. */ +Real thickness = bh; +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec2d(-DL_sponge, -0.25 * DH), Vec2d(DL, 1.25 * DH)); +/** Observation locations for flow field. */ +StdVec observation_locations = {Vec2d(3.0, 5.0), Vec2d(4.0, 5.0), Vec2d(5.0, 5.0)}; +/** + * @brief Material properties of the water and the flexible strucutre. + */ +// Real rho0_f = 1.18e-3; +// Real U_f = 51.3; +// Real c_f = 10.0 * U_f; +// Real Re = 332.6; +// Real mu_f = rho0_f * U_f * (2.0 * insert_circle_radius) / Re; + +// Real rho0_s = 0.1; +// Real poisson = 0.45; +// Real Youngs_modulus =2.5e6; +// Real physical_viscosity = 0.1 * Youngs_modulus; +Real rho0_f = 1.0; +Real U_f = 1.0; +Real c_f = 10.0 * U_f; +Real Re = 100.0; +Real mu_f = rho0_f * U_f * (2.0 * insert_circle_radius) / Re; + +Real rho0_s = 10.0; +Real poisson = 0.45; +Real Ae = 5.0e3; +Real Youngs_modulus = Ae * rho0_f * U_f * U_f; +Real physical_viscosity = 0.1 * Youngs_modulus; +/** Geometry for bodies. */ +std::vector water_block_shape{ Vecd(-DL_sponge, 0.0),Vecd(-DL_sponge, DH), Vecd(DL, DH), Vecd(DL, 0.0), Vecd(-DL_sponge, 0.0)}; +/** create a filament shape */ +Real hbh = bh / 2.0; +Vec2d BLB(insert_circle_center[0], insert_circle_center[1] - hbh); +Vec2d BLT(insert_circle_center[0], insert_circle_center[1] + hbh); +Vec2d BRB(insert_circle_center[0] + insert_circle_radius + bl, insert_circle_center[1] - hbh); +Vec2d BRT(insert_circle_center[0] + insert_circle_radius + bl, insert_circle_center[1] + hbh); +/** create a gate shape */ +std::vector filament_shape{BLB, BLT, BRT, BRB, BLB}; +/** create a square shape */ +Vec2d SLB(insert_circle_center[0]- insert_circle_radius, insert_circle_center[1] - insert_circle_radius); +Vec2d SLT(insert_circle_center[0]- insert_circle_radius, insert_circle_center[1] + insert_circle_radius); +Vec2d SRB(insert_circle_center[0] + insert_circle_radius, insert_circle_center[1] - insert_circle_radius); +Vec2d SRT(insert_circle_center[0] + insert_circle_radius, insert_circle_center[1] + insert_circle_radius); +std::vector square_shape{SLB, SLT, SRT, SRB, SLB}; + +Vec2d emitter_halfsize = Vec2d(0.5 * BW, 0.5 * DH); +Vec2d emitter_translation = Vec2d(-DL_sponge, 0.0) + emitter_halfsize; +Vec2d emitter_buffer_halfsize = Vec2d(0.5 * DL_sponge, 0.5 * DH); +Vec2d emitter_buffer_translation = Vec2d(-DL_sponge, 0.0) + emitter_buffer_halfsize; +Vec2d disposer_halfsize = Vec2d(0.5 * BW, 0.75 * DH); +Vec2d disposer_translation = Vec2d(DL, DH + 0.25 * DH) - disposer_halfsize; +/** Water block shape. */ +class WaterBlock : public MultiPolygonShape +{ +public: + explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(water_block_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(filament_shape, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(square_shape, ShapeBooleanOps::sub); + } +}; +/* Case-dependent geometries.*/ +class WallBoundary : public MultiPolygonShape +{ +public: + explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(square_shape, ShapeBooleanOps::add); + } +}; +/** Particle generator and constraint boundary for shell baffle. */ +int particle_number_mid_surface = int( bl / particle_spacing_ref) + 1; +class ShellBaffleParticleGenerator : public SurfaceParticleGenerator +{ +public: + explicit ShellBaffleParticleGenerator(SPHBody &sph_body) : SurfaceParticleGenerator(sph_body){}; + virtual void initializeGeometricVariables() override + { + for (int i = 0; i < particle_number_mid_surface; i++) + { + Real x = insert_circle_center[0] + insert_circle_radius + Real(i) * particle_spacing_ref - 0.5 * particle_spacing_ref; + Real y = insert_circle_center[1] + 0.5 * particle_spacing_ref; + initializePositionAndVolumetricMeasure(Vecd(x, y), particle_spacing_ref); + Vec2d normal_direction = Vec2d(0, 1.0); + initializeSurfaceProperties(normal_direction, thickness); + } + } +}; +/** Define the boundary geometry. */ +class BoundaryGeometry : public BodyPartByParticle +{ +public: + BoundaryGeometry(SPHBody &body, const std::string &body_part_name) + : BodyPartByParticle(body, body_part_name) + { + TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometry::tagManually, this, _1); + tagParticles(tagging_particle_method); + }; + virtual ~BoundaryGeometry(){}; + +private: + void tagManually(size_t index_i) + { + if (base_particles_.pos_[index_i][0] <= insert_circle_center[0] + insert_circle_radius) + { + body_part_particles_.push_back(index_i); + } + }; +}; +/** + * Free-stream velocity +*/ +struct FreeStreamVelocity +{ + Real u_ref_, t_ref_; + + template + FreeStreamVelocity(BoundaryConditionType &boundary_condition) + : u_ref_(U_f), t_ref_(2.0) {} + + Vecd operator()(Vecd &position, Vecd &velocity) + { + Vecd target_velocity = Vecd::Zero(); + Real run_time = GlobalStaticVariables::physical_time_; + target_velocity[0] = run_time < t_ref_ ? 0.5 * u_ref_ * (1.0 - cos(Pi * run_time / t_ref_)) : u_ref_; + return target_velocity; + } +}; +/** Define time dependent acceleration in x-direction. */ +class TimeDependentAcceleration : public Gravity +{ + Real du_ave_dt_, u_ref_, t_ref_; + +public: + explicit TimeDependentAcceleration(Vecd gravity_vector) + : Gravity(gravity_vector), t_ref_(2.0), u_ref_(U_f), du_ave_dt_(0) {} + + virtual Vecd InducedAcceleration(Vecd &position) override + { + Real run_time_ = GlobalStaticVariables::physical_time_; + du_ave_dt_ = 0.5 * u_ref_ * (Pi / t_ref_) * sin(Pi * run_time_ / t_ref_); + + return run_time_ < t_ref_ ? Vecd(du_ave_dt_, 0.0) : global_acceleration_; + } +}; \ No newline at end of file