diff --git a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.cpp b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.cpp index 7efb8f3013..bbf62889b4 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.cpp @@ -26,13 +26,16 @@ CompressibleFluidInitialCondition::CompressibleFluidInitialCondition(SPHBody &sp E_(*particles_->getVariableDataByName("TotalEnergy")) {} //=================================================================================================// EulerianCompressibleAcousticTimeStepSize:: - EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body) + EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL) : AcousticTimeStepSize(sph_body), rho_(*particles_->getVariableDataByName("Density")), p_(*particles_->getVariableDataByName("Pressure")), vel_(*particles_->getVariableDataByName("Velocity")), smoothing_length_(sph_body.sph_adaptation_->ReferenceSmoothingLength()), - compressible_fluid_(CompressibleFluid(1.0, 1.4)){}; + compressible_fluid_(CompressibleFluid(1.0, 1.4)) +{ + acousticCFL_ = acousticCFL; +}; //=================================================================================================// Real EulerianCompressibleAcousticTimeStepSize::reduce(size_t index_i, Real dt) { @@ -41,7 +44,7 @@ Real EulerianCompressibleAcousticTimeStepSize::reduce(size_t index_i, Real dt) //=================================================================================================// Real EulerianCompressibleAcousticTimeStepSize::outputResult(Real reduced_value) { - return 0.6 / Dimensions * smoothing_length_ / (reduced_value + TinyReal); + return acousticCFL_ / Dimensions * smoothing_length_ / (reduced_value + TinyReal); } //=================================================================================================// } // namespace fluid_dynamics diff --git a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.h b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.h index ea8df6fe47..4426efbfcf 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.h +++ b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_compressible_fluid_integration.h @@ -103,7 +103,7 @@ class EulerianCompressibleAcousticTimeStepSize : public AcousticTimeStepSize Real smoothing_length_; public: - explicit EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body); + explicit EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL = 0.6); virtual ~EulerianCompressibleAcousticTimeStepSize(){}; Real reduce(size_t index_i, Real dt = 0.0); diff --git a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.cpp b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.cpp index 8f12d56400..1b5ed260bf 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.cpp +++ b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.cpp @@ -80,8 +80,26 @@ CompressibleFluidStarState HLLCWithLimiterRiemannSolver:: { Real ul = -e_ij.dot(state_i.vel_); Real ur = -e_ij.dot(state_j.vel_); - Real s_l = ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); - Real s_r = ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); + + //one approach to obtain the smallest and largest wave velocity + //Real s_l = ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); + //Real s_r = ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); + + //another approach to obtain the smallest and largest wave velocity + Vecd vl = state_i.vel_ - ul * (-e_ij); + Vecd vr = state_j.vel_ - ur * (-e_ij); + Real R_lf = state_j.rho_ / state_i.rho_; + Real u_tlide = (ul + ur * R_lf) / (1.0 + R_lf); + Real v_tlide = (vl.norm() + vr.norm() * R_lf) / (1.0 + R_lf); + Real q_tlide = u_tlide; + Real hl = (state_i.E_ + state_i.p_) / state_i.rho_; + Real hr = (state_j.E_ + state_j.p_) / state_j.rho_; + Real h_tlide = (hl + hr * R_lf) / (1.0 + R_lf); + Real sound_tlide = sqrt((1.4 - 1.0) * (h_tlide - 0.5 * (u_tlide * u_tlide + v_tlide * v_tlide))); + Real s_l = SMIN(ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_), q_tlide - sound_tlide); + Real s_r = SMAX(ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_), q_tlide + sound_tlide); + + //obtain the middle wave velocity Real rhol_cl = compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_) * state_i.rho_; Real rhor_cr = compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_) * state_j.rho_; Real clr = (rhol_cl + rhor_cr) / (state_i.rho_ + state_j.rho_); diff --git a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.h b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.h index c87a825540..ab276bc750 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.h +++ b/src/shared/particle_dynamics/fluid_dynamics/eulerian_fluid_dynamics/eulerian_riemann_solver.h @@ -86,7 +86,7 @@ class HLLCWithLimiterRiemannSolver Real limiter_parameter_; public: - HLLCWithLimiterRiemannSolver(CompressibleFluid &compressible_fluid_i, CompressibleFluid &compressible_fluid_j, Real limiter_parameter = 5.0); + HLLCWithLimiterRiemannSolver(CompressibleFluid &compressible_fluid_i, CompressibleFluid &compressible_fluid_j, Real limiter_parameter = 1.0); CompressibleFluidStarState getInterfaceState(const CompressibleFluidState &state_i, const CompressibleFluidState &state_j, const Vecd &e_ij); }; } // namespace SPH diff --git a/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.cpp b/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.cpp new file mode 100644 index 0000000000..ed69a9f697 --- /dev/null +++ b/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.cpp @@ -0,0 +1,162 @@ +#include "eulerian_ghost_boundary.h" + +#include "base_particles.hpp" +namespace SPH +{ +//=================================================================================================// +GhostCreationInESPH::GhostCreationInESPH(BaseInnerRelation &inner_relation, Ghost &ghost_boundary) + : DataDelegateInner(inner_relation), sph_body_(inner_relation.getSPHBody()), + ghost_boundary_(ghost_boundary), get_inner_neighbor_(&inner_relation.getSPHBody()), + indicator_(*particles_->getVariableDataByName("Indicator")), + Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), + pos_(*particles_->getVariableDataByName("Position")), + ghost_bound_(ghost_boundary.GhostBound()) +{ + ghostGenerationAndAddToConfiguration(); +}; +//=================================================================================================// +void GhostCreationInESPH::ghostGenerationAndAddToConfiguration() +{ + ghost_bound_.second = ghost_bound_.first; + RealAndGhostParticleData real_and_ghost_necessary_data; + + for (size_t index_i = 0; index_i != particles_->TotalRealParticles(); ++index_i) + { + if (indicator_[index_i] == 1) + { + Vecd gradient_summation = Vecd::Zero(); + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; + gradient_summation += dW_ijV_j * inner_neighborhood.e_ij_[n]; + } + Real ghost_particle_dW_ijV_j = -gradient_summation.norm(); + Vecd ghost_particle_eij = -gradient_summation / ghost_particle_dW_ijV_j; + Real distance_to_ghost = fabs(sph_body_.getInitialShape().findSignedDistance(pos_[index_i])); + Vecd diaplacement_to_ghost = distance_to_ghost * ghost_particle_eij; + Vecd ghost_position = pos_[index_i] - diaplacement_to_ghost; + mutex_create_ghost_particle_.lock(); + size_t ghost_particle_index = ghost_bound_.second; + ghost_bound_.second++; + ghost_boundary_.checkWithinGhostSize(ghost_bound_); + particles_->updateGhostParticle(ghost_particle_index, index_i); + //Here, we assue the volume as the same the real particle + Vol_[ghost_particle_index] = Vol_[index_i]; + pos_[ghost_particle_index] = ghost_position; + mutex_create_ghost_particle_.unlock(); + Real double_distance_to_ghost = 2.0 * distance_to_ghost; + Real ghost_particle_dW_ij = ghost_particle_dW_ijV_j / (Vol_[ghost_particle_index] + TinyReal); + get_inner_neighbor_(inner_neighborhood, double_distance_to_ghost, ghost_particle_dW_ij, ghost_particle_eij, ghost_particle_index); + + //add all necessary data of real particle and corresponding ghost particle + real_and_ghost_necessary_data.real_index_ = index_i; + real_and_ghost_necessary_data.ghost_index_ = ghost_particle_index; + real_and_ghost_necessary_data.e_ij_ghost_ = ghost_particle_eij; + real_and_ghost_particle_data_.push_back(real_and_ghost_necessary_data); + } + } +} +//=================================================================================================// +GhostBoundaryConditionSetupInESPH:: + GhostBoundaryConditionSetupInESPH(BaseInnerRelation &inner_relation, GhostCreationInESPH &ghost_creation) + : DataDelegateInner(inner_relation), + rho_(*particles_->getVariableDataByName("Density")), + Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), + mass_(*particles_->getVariableDataByName("Mass")), + vel_(*particles_->getVariableDataByName("Velocity")), + pos_(*particles_->getVariableDataByName("Position")), + mom_(*particles_->getVariableDataByName("Momentum")), + ghost_bound_(ghost_creation.ghost_bound_), + real_and_ghost_particle_data_(ghost_creation.real_and_ghost_particle_data_), + boundary_type_(*particles_->registerSharedVariable("BoundaryType")), + W0_(sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd)) +{ + setupBoundaryTypes(); +} +//=================================================================================================// +void GhostBoundaryConditionSetupInESPH::resetBoundaryConditions() +{ + for (size_t ghost_number = 0; ghost_number != real_and_ghost_particle_data_.size(); ++ghost_number) + { + size_t index_i = real_and_ghost_particle_data_[ghost_number].real_index_; + size_t ghost_index = real_and_ghost_particle_data_[ghost_number].ghost_index_; + Vecd e_ig = real_and_ghost_particle_data_[ghost_number].e_ij_ghost_; + + // Dispatch the appropriate boundary condition + switch (boundary_type_[index_i]) + { + case 3: // this refer to the different types of wall boundary conditions + applyNonSlipWallBoundary(ghost_index, index_i); + applyReflectiveWallBoundary(ghost_index, index_i, e_ig); + break; + case 4: + applyTopBoundary(ghost_index, index_i); + break; + case 5: + applyPressureOutletBC(ghost_index, index_i); + break; + case 7: + applySymmetryBoundary(ghost_index, index_i, e_ig); + break; + case 9: + applyFarFieldBoundary(ghost_index, index_i); + break; + case 10: + applyGivenValueInletFlow(ghost_index); + applyVelocityInletFlow(ghost_index, index_i); + break; + case 36: + applyOutletBoundary(ghost_index, index_i); + break; + } + } +} +//=================================================================================================// +GhostKernelGradientUpdate::GhostKernelGradientUpdate(BaseInnerRelation &inner_relation) + : LocalDynamics(inner_relation.getSPHBody()), DataDelegateInner(inner_relation), + Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), + kernel_gradient_original_summation_(*particles_->registerSharedVariable("KernelGradientOriginalSummation")), + indicator_(*particles_->getVariableDataByName("Indicator")){}; +//=================================================================================================// +void GhostKernelGradientUpdate::interaction(size_t index_i, Real dt) +{ + if (indicator_[index_i] == 1) + { + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + Vecd gradient_summation = Vecd::Zero(); + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Vecd &e_ij = inner_neighborhood.e_ij_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; + if (index_j < particles_->TotalRealParticles()) + { + gradient_summation += dW_ijV_j * e_ij; + } + } + kernel_gradient_original_summation_[index_i] = gradient_summation; + } +} +//=================================================================================================// +void GhostKernelGradientUpdate::update(size_t index_i, Real dt) +{ + if (indicator_[index_i] == 1) + { + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + if (index_j >= particles_->TotalRealParticles()) + { + Vecd ghost_dW_ijV_j_with_e_ij = -kernel_gradient_original_summation_[index_i]; + Real ghost_dW_ijV_j = -fabs(ghost_dW_ijV_j_with_e_ij.norm()); + inner_neighborhood.dW_ij_[n] = ghost_dW_ijV_j / Vol_[index_j]; + inner_neighborhood.e_ij_[n] = ghost_dW_ijV_j_with_e_ij / (ghost_dW_ijV_j + TinyReal); + } + } + } +} +//=================================================================================================// +} // namespace SPH diff --git a/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.h b/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.h new file mode 100644 index 0000000000..30d4a6a3b4 --- /dev/null +++ b/src/shared/particle_dynamics/general_dynamics/boundary_condition/eulerian_ghost_boundary.h @@ -0,0 +1,115 @@ +/* ------------------------------------------------------------------------- * + * SPHinXsys * + * ------------------------------------------------------------------------- * + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic * + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4. * + * * + * Portions copyright (c) 2017-2023 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a * + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------- */ +/** + * @file eulerian_ghost_boundary.h + * @brief This is the particle dynamics for domain bounding + * @author Zhentong Wang, Xiangyu Hu + */ +#ifndef Eulerian_GHOST_BOUNDARY_H +#define Eulerian_GHOST_BOUNDARY_H + +#include "base_general_dynamics.h" +#include "particle_reserve.h" +#include "unstructured_mesh.h" +namespace SPH +{ +struct RealAndGhostParticleData +{ + size_t real_index_, ghost_index_; + Vecd e_ij_ghost_; +}; + +/** + * @class GhostCreationInESPH + * @brief The extra storage for the boundary particles, similar with the ghost particle, + * is established to obey the zero-order consistency. + */ +class GhostCreationInESPH : public DataDelegateInner +{ + public: + explicit GhostCreationInESPH(BaseInnerRelation &inner_relation, Ghost &ghost_boundary); + virtual ~GhostCreationInESPH(){}; + std::vector real_and_ghost_particle_data_; + void ghostGenerationAndAddToConfiguration(); + + protected: + SPHBody &sph_body_; + Ghost &ghost_boundary_; + std::mutex mutex_create_ghost_particle_; /**< mutex exclusion for memory conflict */ + NeighborBuilderInnerInFVM get_inner_neighbor_; + StdLargeVec &indicator_; + StdLargeVec &Vol_; + StdLargeVec &pos_; + + public: + std::pair &ghost_bound_; +}; + +//---------------------------------------------------------------------- +// GhostBoundaryConditionSetupInESPH +//---------------------------------------------------------------------- +class GhostBoundaryConditionSetupInESPH : public DataDelegateInner +{ + public: + GhostBoundaryConditionSetupInESPH(BaseInnerRelation &inner_relation, GhostCreationInESPH &ghost_creation); + virtual ~GhostBoundaryConditionSetupInESPH(){}; + virtual void applyReflectiveWallBoundary(size_t ghost_index, size_t index_i, Vecd e_ij){}; + virtual void applyNonSlipWallBoundary(size_t ghost_index, size_t index_i){}; + virtual void applyGivenValueInletFlow(size_t ghost_index){}; + virtual void applyOutletBoundary(size_t ghost_index, size_t index_i){}; + virtual void applyTopBoundary(size_t ghost_index, size_t index_i){}; + virtual void applyFarFieldBoundary(size_t ghost_index, size_t index_i){}; + virtual void applyPressureOutletBC(size_t ghost_index, size_t index_i){}; + virtual void applySymmetryBoundary(size_t ghost_index, size_t index_i, Vecd e_ij){}; + virtual void applyVelocityInletFlow(size_t ghost_index, size_t index_i){}; + virtual void setupBoundaryTypes(){}; + void resetBoundaryConditions(); + + protected: + StdLargeVec &rho_, &Vol_, &mass_; + StdLargeVec &vel_, &pos_, &mom_; + std::pair &ghost_bound_; + std::vector real_and_ghost_particle_data_; + StdLargeVec &boundary_type_; + Real W0_; +}; +/** + * @class GhostKernelGradientUpdate + * @brief Here, this class is used to update the boundary particle + after implementing the kernel correction to strictly follow the zero-order consistency + */ +class GhostKernelGradientUpdate : public LocalDynamics, public DataDelegateInner +{ + public: + explicit GhostKernelGradientUpdate(BaseInnerRelation &inner_relation); + virtual ~GhostKernelGradientUpdate(){}; + void interaction(size_t index_i, Real dt = 0.0); + void update(size_t index_i, Real dt = 0.0); + + protected: + StdLargeVec &Vol_; + StdLargeVec &kernel_gradient_original_summation_; + StdLargeVec &indicator_; +}; +} // namespace SPH +#endif // Eulerian_GHOST_BOUNDARY_H diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.cpp b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.cpp new file mode 100644 index 0000000000..f27260a383 --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.cpp @@ -0,0 +1,183 @@ +/** + * @file 2d_eulerian_supersonic_flow_around_cylinder.cpp + * @brief This is the compressible test for Eulerian supersonic flow around a cylinder + with the FVM boundary algorithm, i.e., zero-order consistency, in the SPHinXsys. + * @author Zhentong Wang and Xiangyu Hu + */ +#include "2d_eulerian_supersonic_flow_around_cylinder.h" +#include "sphinxsys.h" +using namespace SPH; +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem. + //---------------------------------------------------------------------- + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + // Tag for run particle relaxation for the initial body fitted distribution. + sph_system.setRunParticleRelaxation(false); + // Tag for computation start with relaxed body fitted particles distribution. + sph_system.setReloadParticles(false); + // Handle command line arguments and override the tags for particle relaxation and reload. + sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + FluidBody fluid_block(sph_system, makeShared("FluidBlock")); + fluid_block.sph_adaptation_->resetKernel>(20); + fluid_block.defineBodyLevelSetShape(); + fluid_block.defineMaterial(rho_reference, heat_capacity_ratio); + Ghost ghost_boundary(0.5); + (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) + ? fluid_block.generateParticlesWithReserve(ghost_boundary, fluid_block.getName()) + : fluid_block.generateParticlesWithReserve(ghost_boundary); + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + // Note that the same relation should be defined only once. + //---------------------------------------------------------------------- + InnerRelation fluid_block_inner(fluid_block); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (sph_system.RunParticleRelaxation()) + { + //---------------------------------------------------------------------- + // Methods used for particle relaxation. + //---------------------------------------------------------------------- + using namespace relax_dynamics; + SimpleDynamics random_water_body_particles(fluid_block); + BodyStatesRecordingToVtp write_real_body_states(fluid_block); + ReloadParticleIO write_real_body_particle_reload_files(fluid_block); + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(fluid_block_inner); + //---------------------------------------------------------------------- + // Particle relaxation starts here. + //---------------------------------------------------------------------- + random_water_body_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_real_body_states.writeToFile(0); + + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps N = " << ite_p << "\n"; + write_real_body_states.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process finish !" << std::endl; + + write_real_body_particle_reload_files.writeToFile(0); + + return 0; + } + //---------------------------------------------------------------------- + // Define the main numerical methods used in the simulation. + // Note that there may be data dependence on the constructors of these methods. + //---------------------------------------------------------------------- + InteractionWithUpdate>> surface_indicator(fluid_block_inner); + InteractionWithUpdate pressure_relaxation(fluid_block_inner); + InteractionWithUpdate density_and_energy_relaxation(fluid_block_inner); + SimpleDynamics initial_condition(fluid_block); + ReduceDynamics get_fluid_time_step_size(fluid_block, 0.1); + InteractionWithUpdate kernel_correction_matrix(fluid_block_inner); + InteractionDynamics kernel_gradient_update(fluid_block_inner); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + initial_condition.exec(); + surface_indicator.exec(); + //---------------------------------------------------------------------- + // Prepare the ghost particles and its configuration and Boundary conditions setup. + // Ghost kernel gradient update is to make the boundary particles strictly achieve zero-order consistency. + // strictly achieve zero-order consistency + //---------------------------------------------------------------------- + GhostCreationInESPH ghost_creation_for_boundary_condition(fluid_block_inner, ghost_boundary); + SupersonicFlowBoundaryConditionSetup boundary_condition_setup(fluid_block_inner, ghost_creation_for_boundary_condition); + InteractionWithUpdate ghost_kernel_gradient_update(fluid_block_inner); + boundary_condition_setup.resetBoundaryConditions(); + kernel_correction_matrix.exec(); + kernel_gradient_update.exec(); + ghost_kernel_gradient_update.exec(); + //---------------------------------------------------------------------- + // Define the methods for I/O operations and observations of the simulation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp write_real_body_states(sph_system); + RegressionTestEnsembleAverage> write_maximum_speed(fluid_block); + write_real_body_states.addToWrite(fluid_block, "Indicator"); + write_real_body_states.addToWrite(fluid_block, "Velocity"); + write_real_body_states.addToWrite(fluid_block, "Pressure"); + //---------------------------------------------------------------------- + // Setup for time-stepping control + //---------------------------------------------------------------------- + size_t number_of_iterations = 0; + int screen_output_interval = 500; + Real end_time = 40.0; + Real output_interval = 1.0; /**< time stamps for output. */ + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + write_real_body_states.writeToFile(0); + //---------------------------------------------------------------------- + // Main loop starts here. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + while (integration_time < output_interval) + { + Real dt = get_fluid_time_step_size.exec(); + pressure_relaxation.exec(dt); + boundary_condition_setup.resetBoundaryConditions(); + density_and_energy_relaxation.exec(dt); + boundary_condition_setup.resetBoundaryConditions(); + + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + if (number_of_iterations % screen_output_interval == 0) + { + write_maximum_speed.writeToFile(number_of_iterations); + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " dt = " << dt << "\n"; + } + number_of_iterations++; + } + + TickCount t2 = TickCount::now(); + write_real_body_states.writeToFile(); + + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (sph_system.GenerateRegressionData()) + { + write_maximum_speed.generateDataBase(1.0e-3, 1.0e-3); + } + else + { + write_maximum_speed.testResult(); + } + + return 0; +} diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.h b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.h new file mode 100644 index 0000000000..e5f05ffeca --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/2d_eulerian_supersonic_flow_around_cylinder.h @@ -0,0 +1,233 @@ +/** + * @file 2d_eulerian_supersonic_flow_around_cylinder.h + * @brief This is the compressible test for Eulerian supersonic flow around a cylinder + with the FVM boundary algorithm in the SPHinXsys. + * @author Zhentong Wang and Xiangyu Hu + */ +#include "sphinxsys.h" +#include "eulerian_ghost_boundary.h" +using namespace SPH; +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +Real particle_spacing_ref = 1.0 / 7.0; /**< Initial reference particle spacing. */ +Real calculation_circle_radius = 11.0; +Vecd insert_circle_center(7.0, 0.0); +Vecd calculation_circle_center(calculation_circle_radius, 0.0); +Real insert_circle_radius = 1.0; +Real calculation_circle_radius_with_BW = 11.0 + 4.0 * particle_spacing_ref; +BoundingBox system_domain_bounds(Vec2d(0.0, -calculation_circle_radius_with_BW), Vec2d(calculation_circle_radius_with_BW, calculation_circle_radius_with_BW)); +//---------------------------------------------------------------------- +// Material properties of the fluid. +//---------------------------------------------------------------------- +Real rho_reference = 1.0; +Real rho_farfield = 1.0; /**< initial density of one fluid. */ +Real heat_capacity_ratio = 1.4; /**< heat capacity ratio. */ +Real p_farfield = 1.0 / heat_capacity_ratio; /**< initial pressure of one fluid. */ +Real Mach_number = 2.0; +//---------------------------------------------------------------------- +// Define geometries and body shapes +//---------------------------------------------------------------------- +std::vector Creatrightsqureshape() +{ + // geometry + std::vector right_square_shape; + right_square_shape.push_back(Vecd(calculation_circle_center[0], -calculation_circle_radius)); + right_square_shape.push_back(Vecd(calculation_circle_center[0], calculation_circle_radius)); + right_square_shape.push_back(Vecd(calculation_circle_center[0] + calculation_circle_radius, calculation_circle_radius)); + right_square_shape.push_back(Vecd(calculation_circle_center[0] + calculation_circle_radius, -calculation_circle_radius)); + right_square_shape.push_back(Vecd(calculation_circle_center[0], -calculation_circle_radius)); + return right_square_shape; +} +class FluidBlock : public MultiPolygonShape +{ + public: + explicit FluidBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addACircle(calculation_circle_center, calculation_circle_radius, 100, ShapeBooleanOps::add); + multi_polygon_.addACircle(insert_circle_center, insert_circle_radius, 100, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(Creatrightsqureshape(), ShapeBooleanOps::sub); + } +}; +//---------------------------------------------------------------------- +// Case-dependent initial condition. +//---------------------------------------------------------------------- +class SupersonicFlowInitialCondition : public fluid_dynamics::CompressibleFluidInitialCondition +{ + public: + explicit SupersonicFlowInitialCondition(SPHBody &sph_body) + : fluid_dynamics::CompressibleFluidInitialCondition(sph_body){}; + void update(size_t index_i, Real dt) + { + p_[index_i] = p_farfield; + rho_[index_i] = rho_farfield; + mass_[index_i] = rho_[index_i] * Vol_[index_i]; + Real sound_speed = sqrt(p_[index_i] * gamma_ / rho_[index_i]); + vel_[index_i][0] = Mach_number * sound_speed; + vel_[index_i][1] = 0.0; + mom_[index_i] = mass_[index_i] * vel_[index_i]; + Real rho_e = p_[index_i] / (gamma_ - 1.0); + E_[index_i] = rho_e * Vol_[index_i] + 0.5 * mass_[index_i] * vel_[index_i].squaredNorm(); + } + + protected: + Real gamma_ = heat_capacity_ratio; +}; + +//---------------------------------------------------------------------- +// SupersonicFlowBoundaryConditionSetup +//---------------------------------------------------------------------- +class SupersonicFlowBoundaryConditionSetup : public GhostBoundaryConditionSetupInESPH +{ + public: + SupersonicFlowBoundaryConditionSetup(BaseInnerRelation &inner_relation, GhostCreationInESPH &ghost_creation) + : GhostBoundaryConditionSetupInESPH(inner_relation, ghost_creation), + p_(*particles_->getVariableDataByName("Pressure")), + E_(*particles_->getVariableDataByName("TotalEnergy")) + { + setupBoundaryTypes(); + }; + virtual ~SupersonicFlowBoundaryConditionSetup(){}; + + //Here, we follow the FVM labelling different boundary conditions with different number + void setupBoundaryTypes() + { + for (size_t ghost_number = 0; ghost_number != real_and_ghost_particle_data_.size(); ++ghost_number) + { + size_t index_i = real_and_ghost_particle_data_[ghost_number].real_index_; + //around cylinder + if ((pos_[index_i] - insert_circle_center).norm() <= insert_circle_radius + 5.0 * particle_spacing_ref) + { + boundary_type_[index_i] = 3; + } + //outer boundary + else if ((pos_[index_i] - insert_circle_center).norm() > insert_circle_radius + 5.0 * particle_spacing_ref) + { + boundary_type_[index_i] = 9; + } + } + }; + + // For inner boundaries, we set the refective boundary conditions as boundary_type 3. + void applyReflectiveWallBoundary(size_t ghost_index, size_t index_i, Vecd e_ij) override + { + rho_[ghost_index] = rho_[index_i]; + p_[ghost_index] = p_[index_i]; + Real rho_e = p_[ghost_index] / (heat_capacity_ratio - 1.0); + vel_[ghost_index] = (vel_[index_i] - e_ij.dot(vel_[index_i]) * (e_ij)) + (-e_ij.dot(vel_[index_i]) * (e_ij)); + mass_[ghost_index] = rho_[ghost_index] * Vol_[ghost_index]; + mom_[ghost_index] = mass_[ghost_index] * vel_[ghost_index]; + E_[ghost_index] = rho_e * Vol_[ghost_index] + 0.5 * mass_[ghost_index] * vel_[ghost_index].squaredNorm(); + } + + //For outer boundaries, we set the non-refective far-field boundary conditions as boundary_type 9. + void applyFarFieldBoundary(size_t ghost_index, size_t index_i) override + { + Vecd normal_direction_index_i = sph_body_.getInitialShape().findNormalDirection(pos_[index_i]); + Vecd velocity_farfield = Vecd::Zero(); + Real sound_speed = sqrt(p_farfield * heat_capacity_ratio / rho_farfield); + velocity_farfield[0] = Mach_number * sound_speed; + Real velocity_farfield_normal = velocity_farfield.dot(normal_direction_index_i); + Real velocity_boundary_normal = vel_[index_i].dot(normal_direction_index_i); + + //judge the inflow boundary condition + if (normal_direction_index_i[0] <= 0.0 || fabs(normal_direction_index_i[1]) > fabs(normal_direction_index_i[0])) + { + //supersonic inflow condition + if (fabs(velocity_boundary_normal) >= sound_speed) + { + vel_[ghost_index] = velocity_farfield; + p_[ghost_index] = p_farfield; + rho_[ghost_index] = rho_farfield; + mass_[ghost_index] = rho_[ghost_index] * Vol_[ghost_index]; + mom_[ghost_index] = mass_[ghost_index] * vel_[ghost_index]; + Real rho_e = p_[ghost_index] / (heat_capacity_ratio - 1.0); + E_[ghost_index] = rho_e * Vol_[ghost_index] + 0.5 * mass_[ghost_index] * vel_[ghost_index].squaredNorm(); + } + //subsonic inflow condition + if (fabs(velocity_boundary_normal) < sound_speed) + { + Real inner_weight_summation = W0_ * Vol_[index_i]; + Real rho_summation = 0.0; + Real p_summation = 0.0; + Real vel_normal_summation = 0.0; + size_t total_inner_neighbor_particles = 0; + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real W_ij = inner_neighborhood.W_ij_[n]; + inner_weight_summation += W_ij * Vol_[index_j]; + rho_summation += rho_[index_j]; + vel_normal_summation += vel_[index_j].dot(normal_direction_index_i); + p_summation += p_[index_j]; + total_inner_neighbor_particles += 1; + } + Real rho_average = rho_summation / (total_inner_neighbor_particles + TinyReal); + Real vel_normal_average = vel_normal_summation / (total_inner_neighbor_particles + TinyReal); + Real p_average = p_summation / (total_inner_neighbor_particles + TinyReal); + + p_[ghost_index] = p_average * inner_weight_summation + p_farfield * (1.0 - inner_weight_summation); + rho_[ghost_index] = rho_average * inner_weight_summation + rho_farfield * (1.0 - inner_weight_summation); + Real vel_normal = vel_normal_average * inner_weight_summation + velocity_farfield_normal * (1.0 - inner_weight_summation); + vel_[ghost_index] = vel_normal * normal_direction_index_i + (velocity_farfield - velocity_farfield_normal * normal_direction_index_i); + mass_[ghost_index] = rho_[ghost_index] * Vol_[ghost_index]; + mom_[ghost_index] = mass_[ghost_index] * vel_[ghost_index]; + Real rho_e = p_[ghost_index] / (heat_capacity_ratio - 1.0); + E_[ghost_index] = rho_e * Vol_[ghost_index] + 0.5 * mass_[ghost_index] * vel_[ghost_index].squaredNorm(); + } + } + else + { + //supersonic outflow condition + if (fabs(velocity_boundary_normal) >= sound_speed) + { + vel_[ghost_index] = vel_[index_i]; + p_[ghost_index] = p_[index_i]; + rho_[ghost_index] = rho_[index_i]; + mass_[ghost_index] = rho_[ghost_index] * Vol_[ghost_index]; + mom_[ghost_index] = mass_[ghost_index] * vel_[ghost_index]; + Real rho_e = p_[ghost_index] / (heat_capacity_ratio - 1.0); + E_[ghost_index] = rho_e * Vol_[ghost_index] + 0.5 * mass_[ghost_index] * vel_[ghost_index].squaredNorm(); + } + //subsonic outflow condition + if (fabs(velocity_boundary_normal) < sound_speed) + { + Real inner_weight_summation = W0_ * Vol_[index_i]; + Real rho_summation = 0.0; + Real p_summation = 0.0; + Real vel_normal_summation(0.0); + Vecd vel_tangential_summation = Vecd::Zero(); + size_t total_inner_neighbor_particles = 0; + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real W_ij = inner_neighborhood.W_ij_[n]; + inner_weight_summation += W_ij * Vol_[index_j]; + rho_summation += rho_[index_j]; + vel_normal_summation += vel_[index_j].dot(normal_direction_index_i); + vel_tangential_summation += vel_[index_j] - vel_[index_j].dot(normal_direction_index_i) * normal_direction_index_i; + p_summation += p_[index_j]; + total_inner_neighbor_particles += 1; + } + Real rho_average = rho_summation / (total_inner_neighbor_particles + TinyReal); + Real vel_normal_average = vel_normal_summation / (total_inner_neighbor_particles + TinyReal); + Vecd vel_tangential_average = vel_tangential_summation / (total_inner_neighbor_particles + TinyReal); + Real p_average = p_summation / (total_inner_neighbor_particles + TinyReal); + + p_[ghost_index] = p_average * inner_weight_summation + p_farfield * (1.0 - inner_weight_summation); + rho_[ghost_index] = rho_average * inner_weight_summation + rho_farfield * (1.0 - inner_weight_summation); + Real vel_normal = vel_normal_average * inner_weight_summation + velocity_farfield_normal * (1.0 - inner_weight_summation); + vel_[ghost_index] = vel_normal * normal_direction_index_i + vel_tangential_average; + mass_[ghost_index] = rho_[ghost_index] * Vol_[ghost_index]; + mom_[ghost_index] = mass_[ghost_index] * vel_[ghost_index]; + Real rho_e = p_[ghost_index] / (heat_capacity_ratio - 1.0); + E_[ghost_index] = rho_e * Vol_[ghost_index] + 0.5 * mass_[ghost_index] * vel_[ghost_index].squaredNorm(); + } + } + } + + protected: + StdLargeVec &p_, &E_; +}; diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/CMakeLists.txt b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/CMakeLists.txt new file mode 100644 index 0000000000..21c4416b22 --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/CMakeLists.txt @@ -0,0 +1,24 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ + DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +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}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "periodic boundary, Eulerian") diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_ensemble_averaged_mean_variance.xml b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_ensemble_averaged_mean_variance.xml new file mode 100644 index 0000000000..45c43020c7 --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_ensemble_averaged_mean_variance.xml @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_runtimes.dat b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_runtimes.dat new file mode 100644 index 0000000000..fd93545d1b --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/FluidBlock_MaximumSpeed_runtimes.dat @@ -0,0 +1,3 @@ +true +8 +5 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/regression_test_tool.py b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..8d1f22017c --- /dev/null +++ b/tests/2d_examples/test_2d_eulerian_supersonic_flow_new_BC/regression_test_tool/regression_test_tool.py @@ -0,0 +1,36 @@ +# !/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: 2d_eulerian_supersonic_flow_around_cylinder +""" + +case_name = "2d_eulerian_supersonic_flow_around_cylinder" +body_name = "FluidBlock" +parameter_name = "MaximumSpeed" + +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