diff --git a/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h index a4d246e0bd..ceff6d0f03 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h +++ b/src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h @@ -120,7 +120,7 @@ using TransportVelocityCorrectionComplex = BaseTransportVelocityCorrectionComplex; template -using TransportVelocityCorrectionCorrectedComplex = +using TransportVelocityCorrectionCorrectedComplex = BaseTransportVelocityCorrectionComplex; template @@ -128,10 +128,9 @@ using TransportVelocityLimitedCorrectionComplex = BaseTransportVelocityCorrectionComplex; template -using TransportVelocityKimitedCorrectionCorrectedComplex = +using TransportVelocityLimitedCorrectionCorrectedComplex = BaseTransportVelocityCorrectionComplex; - template using TransportVelocityCorrectionComplexAdaptive = BaseTransportVelocityCorrectionComplex; diff --git a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.h b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.h index 2a2b47df9c..099c4a75b3 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.h +++ b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.h @@ -80,8 +80,8 @@ class ViscousForce Real smoothing_length_; }; -template -class ViscousForce, ViscosityType> +template +class ViscousForce, ViscosityType, KernelCorrectionType> : public ViscousForce, public ForcePrior { public: @@ -91,43 +91,43 @@ class ViscousForce, ViscosityType> protected: ViscosityType mu_; + KernelCorrectionType kernel_correction_; }; -using ViscousForceInner = ViscousForce, FixedViscosity>; +using ViscousForceInner = ViscousForce, FixedViscosity, NoKernelCorrection>; -template -class ViscousForce, ViscosityType> +template +class ViscousForce, ViscosityType, KernelCorrectionType> : public ViscousForce, public ForcePrior { public: - explicit ViscousForce(BaseInnerRelation &inner_relation) - : ViscousForce(inner_relation), - ForcePrior(particles_, "ViscousForce"), - mu_(particles_){}; + explicit ViscousForce(BaseInnerRelation &inner_relation); virtual ~ViscousForce(){}; void interaction(size_t index_i, Real dt = 0.0); protected: ViscosityType mu_; + KernelCorrectionType kernel_correction_; }; using BaseViscousForceWithWall = InteractionWithWall; -template -class ViscousForce, ViscosityType> : public BaseViscousForceWithWall +template +class ViscousForce, ViscosityType, KernelCorrectionType> + : public BaseViscousForceWithWall { public: - explicit ViscousForce(BaseContactRelation &wall_contact_relation) - : BaseViscousForceWithWall(wall_contact_relation), - mu_(particles_){}; + explicit ViscousForce(BaseContactRelation &wall_contact_relation); virtual ~ViscousForce(){}; void interaction(size_t index_i, Real dt = 0.0); protected: ViscosityType mu_; + KernelCorrectionType kernel_correction_; }; -template -class ViscousForce, ViscosityType> : public BaseViscousForceWithWall +template +class ViscousForce, ViscosityType, KernelCorrectionType> + : public BaseViscousForceWithWall { public: explicit ViscousForce(BaseContactRelation &wall_contact_relation); @@ -136,11 +136,12 @@ class ViscousForce, ViscosityType> : public B protected: ViscosityType mu_; - StdLargeVec &distance_from_wall_; + KernelCorrectionType kernel_correction_; }; -template -class ViscousForce, ViscosityType> : public ViscousForce +template +class ViscousForce, ViscosityType, KernelCorrectionType> + : public ViscousForce { public: explicit ViscousForce(BaseContactRelation &contact_relation); @@ -149,16 +150,19 @@ class ViscousForce, ViscosityType> : public ViscousForce contact_mu_; + KernelCorrectionType kernel_correction_; + StdVec contact_kernel_corrections_; StdVec *> contact_vel_; StdVec *> wall_Vol_; }; -using ViscousForceWithWall = ComplexInteraction, Contact>, FixedViscosity>; -using MultiPhaseViscousForceWithWall = ComplexInteraction, Contact<>, Contact>, FixedViscosity>; +using ViscousForceWithWall = ComplexInteraction, Contact>, FixedViscosity, NoKernelCorrection>; +using ViscousForceWithWallCorrection = ComplexInteraction, Contact>, FixedViscosity, LinearGradientCorrection>; +using MultiPhaseViscousForceWithWall = ComplexInteraction, Contact<>, Contact>, FixedViscosity, NoKernelCorrection>; template using NonNewtonianViscousForceWithWall = - ComplexInteraction, Contact>, VariableViscosity>; + ComplexInteraction, Contact>, VariableViscosity, NoKernelCorrection>; /** * @class VorticityInner diff --git a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp index acbe3caaca..cc65c8f7e5 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp @@ -18,17 +18,19 @@ ViscousForce::ViscousForce(BaseRelationType &base_relation) viscous_force_(*this->particles_->template registerSharedVariable("ViscousForce")), smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength()) {} //=================================================================================================// -template -ViscousForce, ViscosityType>::ViscousForce(BaseInnerRelation &inner_relation) +template +ViscousForce, ViscosityType, KernelCorrectionType>:: + ViscousForce(BaseInnerRelation &inner_relation) : ViscousForce(inner_relation), - ForcePrior(particles_, "ViscousForce"), mu_(particles_) + ForcePrior(particles_, "ViscousForce"), mu_(particles_), + kernel_correction_(particles_) { static_assert(std::is_base_of::value, "ParticleAverage is not the base of ViscosityType!"); } //=================================================================================================// -template -void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) +template +void ViscousForce, ViscosityType, KernelCorrectionType>::interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; @@ -39,15 +41,23 @@ void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) // viscous force Vecd vel_derivative = (vel_[index_i] - vel_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); - force += 2.0 * mu_(index_i, index_j) * vel_derivative * + force += (kernel_correction_(index_i) + kernel_correction_(index_j)) * + mu_(index_i, index_j) * vel_derivative * inner_neighborhood.dW_ij_[n] * Vol_[index_j]; } viscous_force_[index_i] = force * Vol_[index_i]; } //=================================================================================================// -template -void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) +template +ViscousForce, ViscosityType, KernelCorrectionType>:: + ViscousForce(BaseInnerRelation &inner_relation) + : ViscousForce(inner_relation), ForcePrior(particles_, "ViscousForce"), + mu_(particles_), kernel_correction_(particles_) {} +//=================================================================================================// +template +void ViscousForce, ViscosityType, KernelCorrectionType>:: + interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; @@ -60,16 +70,24 @@ void ViscousForce, ViscosityType>::interaction(size_t /** The following viscous force is given in Monaghan 2005 (Rep. Prog. Phys.), it seems that * this formulation is more accurate than the previous one for Taylor-Green-Vortex flow. */ Real v_r_ij = (vel_[index_i] - vel_[index_j]).dot(e_ij); - Real eta_ij = 2.0 * Real(Dimensions + 2) * mu_(index_i, index_j) * v_r_ij / + Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_j) * v_r_ij / (r_ij + 0.01 * smoothing_length_); - force += eta_ij * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * e_ij; + force += (kernel_correction_(index_i) + kernel_correction_(index_j)) * + eta_ij * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * e_ij; } viscous_force_[index_i] = force * Vol_[index_i]; } //=================================================================================================// -template -void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) +template +ViscousForce, ViscosityType, KernelCorrectionType>:: + ViscousForce(BaseContactRelation &wall_contact_relation) + : BaseViscousForceWithWall(wall_contact_relation), + mu_(particles_), kernel_correction_(particles_) {} +//=================================================================================================// +template +void ViscousForce, ViscosityType, KernelCorrectionType>:: + interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); for (size_t k = 0; k < contact_configuration_.size(); ++k) @@ -84,26 +102,25 @@ void ViscousForce, ViscosityType>::interaction(size_t index_i, Rea Vecd vel_derivative = 2.0 * (vel_[index_i] - vel_ave_k[index_j]) / (r_ij + 0.01 * smoothing_length_); - force += 2.0 * mu_(index_i, index_i) * vel_derivative * - contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j]; + force += 2.0 * kernel_correction_(index_i) * mu_(index_i, index_i) * + vel_derivative * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j]; } } viscous_force_[index_i] += force * Vol_[index_i]; } //=================================================================================================// -template -ViscousForce, ViscosityType>:: +template +ViscousForce, ViscosityType, KernelCorrectionType>:: ViscousForce(BaseContactRelation &wall_contact_relation) : BaseViscousForceWithWall(wall_contact_relation), - mu_(particles_), - distance_from_wall_(*this->particles_->template registerSharedVariable("DistanceFromWall")){}; + mu_(particles_), kernel_correction_(particles_){}; //=================================================================================================// -template -void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) +template +void ViscousForce, ViscosityType, KernelCorrectionType>:: + interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); - const Vecd &distance_from_wall = distance_from_wall_[index_i]; for (size_t k = 0; k < contact_configuration_.size(); ++k) { StdLargeVec &vel_ave_k = *(wall_vel_ave_[k]); @@ -115,32 +132,35 @@ void ViscousForce, ViscosityType>::interactio Vecd &e_ij = contact_neighborhood.e_ij_[n]; Real r_ij = contact_neighborhood.r_ij_[n]; - Vecd distance_diff = distance_from_wall - r_ij * e_ij; - Real factor = 1.0 - distance_from_wall.dot(distance_diff) / distance_from_wall.squaredNorm(); - Real v_r_ij = factor * (vel_[index_i] - vel_ave_k[index_j]).dot(e_ij); - Real eta_ij = 2.0 * Real(Dimensions + 2) * mu_(index_i, index_i) * v_r_ij / + Real v_r_ij = 2.0 * (vel_[index_i] - vel_ave_k[index_j]).dot(e_ij); + Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_i) * v_r_ij / (r_ij + 0.01 * smoothing_length_); - force += eta_ij * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij; + force += 2.0 * kernel_correction_(index_i) * eta_ij * + contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij; } } viscous_force_[index_i] += force * Vol_[index_i]; } //=================================================================================================// -template -ViscousForce, ViscosityType>::ViscousForce(BaseContactRelation &contact_relation) - : ViscousForce(contact_relation) +template +ViscousForce, ViscosityType, KernelCorrectionType>:: + ViscousForce(BaseContactRelation &contact_relation) + : ViscousForce(contact_relation), + kernel_correction_(particles_) { for (size_t k = 0; k != contact_particles_.size(); ++k) { contact_mu_.emplace_back(ViscosityType(particles_, contact_particles_[k])); + contact_kernel_corrections_.emplace_back(KernelCorrectionType(contact_particles_[k])); contact_vel_.push_back(contact_particles_[k]->template getVariableDataByName("Velocity")); wall_Vol_.push_back(contact_particles_[k]->template getVariableDataByName("VolumetricMeasure")); } } //=================================================================================================// -template -void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt) +template +void ViscousForce, ViscosityType, KernelCorrectionType>:: + interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); for (size_t k = 0; k < contact_configuration_.size(); ++k) @@ -148,13 +168,15 @@ void ViscousForce, ViscosityType>::interaction(size_t index_i, Real dt auto &contact_mu_k = contact_mu_[k]; StdLargeVec &vel_k = *(contact_vel_[k]); StdLargeVec &wall_Vol_k = *(wall_Vol_[k]); + KernelCorrectionType &kernel_correction_k = contact_kernel_corrections_[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 = (vel_[index_i] - vel_k[index_j]) / (contact_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); - force += 2.0 * contact_mu_k(index_i, index_j) * vel_derivative * + force += (kernel_correction_(index_i) + kernel_correction_k(index_j)) * + contact_mu_k(index_i, index_j) * vel_derivative * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j]; } } diff --git a/tests/2d_examples/test_2d_lid_driven_cavity_corrected/lid_driven_cavity.cpp b/tests/2d_examples/test_2d_lid_driven_cavity_corrected/lid_driven_cavity.cpp index 4fb3243174..d71dba4172 100644 --- a/tests/2d_examples/test_2d_lid_driven_cavity_corrected/lid_driven_cavity.cpp +++ b/tests/2d_examples/test_2d_lid_driven_cavity_corrected/lid_driven_cavity.cpp @@ -9,276 +9,277 @@ using namespace SPH; // Namespace cite here. //---------------------------------------------------------------------- // Basic geometry parameters and numerical setup. //---------------------------------------------------------------------- -Real DL = 1.0; /**< box length. */ -Real DH = 1.0; /**< box height. */ +Real DL = 1.0; /**< box length. */ +Real DH = 1.0; /**< box height. */ Real resolution_ref = 1.0 / 50.0; /**< Global reference resolution. */ -Real BW = resolution_ref * 6; /**< Extending width for BCs. */ +Real BW = resolution_ref * 6; /**< Extending width for BCs. */ /** Domain bounds of the system. */ BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(DL + BW, DH + BW)); //---------------------------------------------------------------------- // Material properties of the fluid. //---------------------------------------------------------------------- -Real rho0_f = 1.0; /**< Reference density of fluid. */ -Real U_f = 1.0; /**< Characteristic velocity. */ -Real c_f = 10.0 * U_f; /**< Reference sound speed. */ -Real Re = 100.0; /**< Reynolds number. */ +Real rho0_f = 1.0; /**< Reference density of fluid. */ +Real U_f = 1.0; /**< Characteristic velocity. */ +Real c_f = 10.0 * U_f; /**< Reference sound speed. */ +Real Re = 100.0; /**< Reynolds number. */ Real mu_f = rho0_f * U_f * DL / Re; /**< Dynamics viscosity. */ //---------------------------------------------------------------------- // Cases-dependent geometries //---------------------------------------------------------------------- class WaterBlock : public MultiPolygonShape { -public: - explicit WaterBlock(const std::string& shape_name) : MultiPolygonShape(shape_name) - { - /** Geometry definition. */ - std::vector water_body_shape; - water_body_shape.push_back(Vecd(0.0, 0.0)); - water_body_shape.push_back(Vecd(0.0, DH)); - water_body_shape.push_back(Vecd(DL, DH)); - water_body_shape.push_back(Vecd(DL, 0.0)); - water_body_shape.push_back(Vecd(0.0, 0.0)); - multi_polygon_.addAPolygon(water_body_shape, ShapeBooleanOps::add); - } + public: + explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + /** Geometry definition. */ + std::vector water_body_shape; + water_body_shape.push_back(Vecd(0.0, 0.0)); + water_body_shape.push_back(Vecd(0.0, DH)); + water_body_shape.push_back(Vecd(DL, DH)); + water_body_shape.push_back(Vecd(DL, 0.0)); + water_body_shape.push_back(Vecd(0.0, 0.0)); + multi_polygon_.addAPolygon(water_body_shape, ShapeBooleanOps::add); + } }; /** * @brief Wall boundary body definition. */ class WallBoundary : public MultiPolygonShape { -public: - explicit WallBoundary(const std::string& shape_name) : MultiPolygonShape(shape_name) - { - /** Geometry definition. */ - std::vector outer_wall_shape; - outer_wall_shape.push_back(Vecd(-BW, -BW)); - outer_wall_shape.push_back(Vecd(-BW, DH + BW)); - outer_wall_shape.push_back(Vecd(DL + BW, DH + BW)); - outer_wall_shape.push_back(Vecd(DL + BW, -BW)); - outer_wall_shape.push_back(Vecd(-BW, -BW)); - std::vector inner_wall_shape; - inner_wall_shape.push_back(Vecd(0.0, 0.0)); - inner_wall_shape.push_back(Vecd(0.0, DH)); - inner_wall_shape.push_back(Vecd(DL, DH)); - inner_wall_shape.push_back(Vecd(DL, 0.0)); - inner_wall_shape.push_back(Vecd(0.0, 0.0)); + public: + explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + /** Geometry definition. */ + std::vector outer_wall_shape; + outer_wall_shape.push_back(Vecd(-BW, -BW)); + outer_wall_shape.push_back(Vecd(-BW, DH + BW)); + outer_wall_shape.push_back(Vecd(DL + BW, DH + BW)); + outer_wall_shape.push_back(Vecd(DL + BW, -BW)); + outer_wall_shape.push_back(Vecd(-BW, -BW)); + std::vector inner_wall_shape; + inner_wall_shape.push_back(Vecd(0.0, 0.0)); + inner_wall_shape.push_back(Vecd(0.0, DH)); + inner_wall_shape.push_back(Vecd(DL, DH)); + inner_wall_shape.push_back(Vecd(DL, 0.0)); + inner_wall_shape.push_back(Vecd(0.0, 0.0)); - multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); - } + multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); + } }; //---------------------------------------------------------------------- // Application dependent initial condition //---------------------------------------------------------------------- class BoundaryVelocity : public MotionConstraint { -public: - BoundaryVelocity(SPHBody& body) - : MotionConstraint(body) {} + public: + BoundaryVelocity(SPHBody &body) + : MotionConstraint(body) {} - void update(size_t index_i, Real dt = 0.0) - { - if (pos_[index_i][1] > DH) - { - vel_[index_i][0] = 1.0; - vel_[index_i][1] = 0.0; - } - }; + void update(size_t index_i, Real dt = 0.0) + { + if (pos_[index_i][1] > DH) + { + vel_[index_i][0] = 1.0; + vel_[index_i][1] = 0.0; + } + }; }; //---------------------------------------------------------------------- // An observer particle generator. //---------------------------------------------------------------------- StdVec VelocityXObserverParticle() { - StdVec observation_points; - size_t number_of_observation_point = 5; - Real range_of_measure = 1.0 - 0.5 * resolution_ref; - Real start_of_measure = 0.5 * resolution_ref; + StdVec observation_points; + size_t number_of_observation_point = 5; + Real range_of_measure = 1.0 - 0.5 * resolution_ref; + Real start_of_measure = 0.5 * resolution_ref; - for (size_t i = 0; i < number_of_observation_point; ++i) - { - Vec2d point_corrdinate(range_of_measure * (Real)i / (Real)(number_of_observation_point - 1) + start_of_measure, 0.5 * DL); - observation_points.push_back(point_corrdinate); - } - return observation_points; + for (size_t i = 0; i < number_of_observation_point; ++i) + { + Vec2d point_corrdinate(range_of_measure * (Real)i / (Real)(number_of_observation_point - 1) + start_of_measure, 0.5 * DL); + observation_points.push_back(point_corrdinate); + } + return observation_points; } StdVec VelocityYObserverParticle() { - StdVec observation_points; - size_t number_of_observation_point = 5; - Real range_of_measure = 1.0 - 0.5 * resolution_ref; - Real start_of_measure = 0.5 * resolution_ref; - for (size_t i = 0; i < number_of_observation_point; ++i) - { - Vec2d point_corrdinate(0.5 * DH, range_of_measure * (Real)i / - (Real)(number_of_observation_point - 1) + start_of_measure); - observation_points.push_back(point_corrdinate); - } - return observation_points; + StdVec observation_points; + size_t number_of_observation_point = 5; + Real range_of_measure = 1.0 - 0.5 * resolution_ref; + Real start_of_measure = 0.5 * resolution_ref; + for (size_t i = 0; i < number_of_observation_point; ++i) + { + Vec2d point_corrdinate(0.5 * DH, range_of_measure * (Real)i / + (Real)(number_of_observation_point - 1) + + start_of_measure); + observation_points.push_back(point_corrdinate); + } + return observation_points; } //---------------------------------------------------------------------- // Main program starts here. //---------------------------------------------------------------------- int main(int ac, char *av[]) { - //---------------------------------------------------------------------- - // Build up the environment of a SPHSystem. - //---------------------------------------------------------------------- - SPHSystem sph_system(system_domain_bounds, resolution_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); - /** Set the starting time. */ - GlobalStaticVariables::physical_time_ = 0.0; - IOEnvironment io_environment(sph_system); - sph_system.handleCommandlineOptions(ac, av); - //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - FluidBody water_body(sph_system, makeShared("WaterBody")); - water_body.defineMaterial(rho0_f, c_f, mu_f); - water_body.generateParticles(); + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem. + //---------------------------------------------------------------------- + SPHSystem sph_system(system_domain_bounds, resolution_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); + /** Set the starting time. */ + GlobalStaticVariables::physical_time_ = 0.0; + IOEnvironment io_environment(sph_system); + sph_system.handleCommandlineOptions(ac, av); + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + FluidBody water_body(sph_system, makeShared("WaterBody")); + water_body.defineMaterial(rho0_f, c_f, mu_f); + water_body.generateParticles(); - SolidBody wall_boundary(sph_system, makeShared("Wall")); - wall_boundary.defineMaterial(); - wall_boundary.generateParticles(); - //---------------------------------------------------------------------- - // Particle and body creation of fluid observers. - //---------------------------------------------------------------------- - ObserverBody horizontal_observer(sph_system, "HorizontalVelocity"); - horizontal_observer.generateParticles(VelocityXObserverParticle()); - ObserverBody vertical_observer(sph_system, "VerticalVelocity"); - vertical_observer.generateParticles(VelocityYObserverParticle()); - //---------------------------------------------------------------------- - // 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. - //---------------------------------------------------------------------- - InnerRelation water_block_inner(water_body); - ContactRelation water_block_contact(water_body, { &wall_boundary }); - ContactRelation horizontal_observer_contact(horizontal_observer, { &water_body }); - ContactRelation vertical_observer_contact(vertical_observer, { &water_body }); - ComplexRelation water_block_complex(water_block_inner, water_block_contact); - //---------------------------------------------------------------------- - // Define the main numerical methods used in the simulation. - // Note that there may be data dependence on the constructors of these methods. - //---------------------------------------------------------------------- - SimpleDynamics wall_boundary_normal_direction(wall_boundary); - /** Initial condition with momentum field */ - SimpleDynamics solid_initial_condition(wall_boundary); - /** Kernel correction matrix and transport velocity formulation. */ - InteractionWithUpdate kernel_correction_complex(water_block_inner, water_block_contact); - /** Evaluation of density by summation approach. */ - InteractionWithUpdate update_density_by_summation(water_block_inner, water_block_contact); - /** Pressure and density relaxation algorithm by using verlet time stepping. */ - Dynamics1Level pressure_relaxation(water_block_inner, water_block_contact); - Dynamics1Level density_relaxation(water_block_inner, water_block_contact); - InteractionWithUpdate> - transport_velocity_correction(water_block_inner, water_block_contact); - /** Time step size with considering sound wave speed. */ - ReduceDynamics get_fluid_advection_time_step_size(water_body, U_f); - ReduceDynamics get_fluid_time_step_size(water_body); - /** Computing viscous acceleration with wall. */ - InteractionWithUpdate viscous_acceleration(water_block_inner, water_block_contact); - //---------------------------------------------------------------------- - // Define the methods for I/O operations and observations of the simulation. - //---------------------------------------------------------------------- - /** Output the body states. */ - BodyStatesRecordingToVtp write_real_body_states(sph_system); - write_real_body_states.addToWrite(water_body, "Velocity"); - RegressionTestDynamicTimeWarping> write_horizontal_velocity("Velocity", horizontal_observer_contact); - RegressionTestDynamicTimeWarping> write_vertical_velocity("Velocity", vertical_observer_contact); - //---------------------------------------------------------------------- - // Prepare the simulation with cell linked list, configuration - // and case specified initial condition if necessary. - //---------------------------------------------------------------------- - sph_system.initializeSystemCellLinkedLists(); - sph_system.initializeSystemConfigurations(); - solid_initial_condition.exec(); - write_real_body_states.writeToFile(); - kernel_correction_complex.exec(); - //---------------------------------------------------------------------- - // Setup for time-stepping control - //---------------------------------------------------------------------- - size_t number_of_iterations = 0; - int screen_output_interval = 100; - Real End_Time = 30.0; /**< End time. */ - Real output_interval = 1.0; - Real dt = 1.0; /**< Time stamps for output of body states. */ - //---------------------------------------------------------------------- - // Statistics for CPU time - //---------------------------------------------------------------------- - TickCount t1 = TickCount::now(); + SolidBody wall_boundary(sph_system, makeShared("Wall")); + wall_boundary.defineMaterial(); + wall_boundary.generateParticles(); + //---------------------------------------------------------------------- + // Particle and body creation of fluid observers. + //---------------------------------------------------------------------- + ObserverBody horizontal_observer(sph_system, "HorizontalVelocity"); + horizontal_observer.generateParticles(VelocityXObserverParticle()); + ObserverBody vertical_observer(sph_system, "VerticalVelocity"); + vertical_observer.generateParticles(VelocityYObserverParticle()); + //---------------------------------------------------------------------- + // 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. + //---------------------------------------------------------------------- + InnerRelation water_block_inner(water_body); + ContactRelation water_block_contact(water_body, {&wall_boundary}); + ContactRelation horizontal_observer_contact(horizontal_observer, {&water_body}); + ContactRelation vertical_observer_contact(vertical_observer, {&water_body}); + ComplexRelation water_block_complex(water_block_inner, water_block_contact); + //---------------------------------------------------------------------- + // Define the main numerical methods used in the simulation. + // Note that there may be data dependence on the constructors of these methods. + //---------------------------------------------------------------------- + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + /** Initial condition with momentum field */ + SimpleDynamics solid_initial_condition(wall_boundary); + /** Kernel correction matrix and transport velocity formulation. */ + InteractionWithUpdate kernel_correction_complex(water_block_inner, water_block_contact); + /** Evaluation of density by summation approach. */ + InteractionWithUpdate update_density_by_summation(water_block_inner, water_block_contact); + /** Pressure and density relaxation algorithm by using verlet time stepping. */ + Dynamics1Level pressure_relaxation(water_block_inner, water_block_contact); + Dynamics1Level density_relaxation(water_block_inner, water_block_contact); + InteractionWithUpdate> + transport_velocity_correction(water_block_inner, water_block_contact); + /** Time step size with considering sound wave speed. */ + ReduceDynamics get_fluid_advection_time_step_size(water_body, U_f); + ReduceDynamics get_fluid_time_step_size(water_body); + /** Computing viscous acceleration with wall. */ + InteractionWithUpdate viscous_acceleration(water_block_inner, water_block_contact); + //---------------------------------------------------------------------- + // Define the methods for I/O operations and observations of the simulation. + //---------------------------------------------------------------------- + /** Output the body states. */ + BodyStatesRecordingToVtp write_real_body_states(sph_system); + write_real_body_states.addToWrite(water_body, "Velocity"); + RegressionTestDynamicTimeWarping> write_horizontal_velocity("Velocity", horizontal_observer_contact); + RegressionTestDynamicTimeWarping> write_vertical_velocity("Velocity", vertical_observer_contact); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + solid_initial_condition.exec(); + write_real_body_states.writeToFile(); + kernel_correction_complex.exec(); + //---------------------------------------------------------------------- + // Setup for time-stepping control + //---------------------------------------------------------------------- + size_t number_of_iterations = 0; + int screen_output_interval = 100; + Real End_Time = 30.0; /**< End time. */ + Real output_interval = 1.0; + Real dt = 1.0; /**< Time stamps for output of body states. */ + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); TimeInterval interval; - //---------------------------------------------------------------------- - // First output before the main loop. - //---------------------------------------------------------------------- - /** Output the start states of bodies. */ - 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_advection_time_step_size.exec(); - update_density_by_summation.exec(); - viscous_acceleration.exec(); - - kernel_correction_complex.exec(); - transport_velocity_correction.exec(); - - Real relaxation_time = 0.0; - while(relaxation_time < Dt) - { - // avoid possible smaller acoustic time step size for viscous flow - dt = SMIN(get_fluid_time_step_size.exec(), Dt); - relaxation_time += dt; - integration_time += dt; - pressure_relaxation.exec(dt); - density_relaxation.exec(dt); - GlobalStaticVariables::physical_time_ += dt; - } - - 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 << "\n"; - } - number_of_iterations++; - water_body.updateCellLinkedList(); - water_block_complex.updateConfiguration(); - } - TickCount t2 = TickCount::now(); - write_real_body_states.writeToFile(); - horizontal_observer_contact.updateConfiguration(); - vertical_observer_contact.updateConfiguration(); - write_horizontal_velocity.writeToFile(number_of_iterations); - write_vertical_velocity.writeToFile(number_of_iterations); - TickCount t3 = TickCount::now(); + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + /** Output the start states of bodies. */ + 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_advection_time_step_size.exec(); + update_density_by_summation.exec(); + viscous_acceleration.exec(); + + kernel_correction_complex.exec(); + transport_velocity_correction.exec(); + + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + // avoid possible smaller acoustic time step size for viscous flow + dt = SMIN(get_fluid_time_step_size.exec(), Dt); + relaxation_time += dt; + integration_time += dt; + pressure_relaxation.exec(dt); + density_relaxation.exec(dt); + GlobalStaticVariables::physical_time_ += dt; + } + + 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 << "\n"; + } + number_of_iterations++; + water_body.updateCellLinkedList(); + water_block_complex.updateConfiguration(); + } + TickCount t2 = TickCount::now(); + write_real_body_states.writeToFile(); + horizontal_observer_contact.updateConfiguration(); + vertical_observer_contact.updateConfiguration(); + write_horizontal_velocity.writeToFile(number_of_iterations); + write_vertical_velocity.writeToFile(number_of_iterations); + TickCount t3 = TickCount::now(); interval += t3 - t2; - } - - TickCount t4 = TickCount::now(); + } + + 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_horizontal_velocity.generateDataBase(1.0e-3); - write_vertical_velocity.generateDataBase(1.0e-3); - } - else if (sph_system.RestartStep() == 0) - { - write_horizontal_velocity.testResult(); - write_vertical_velocity.testResult(); - } + if (sph_system.GenerateRegressionData()) + { + write_horizontal_velocity.generateDataBase(1.0e-3); + write_vertical_velocity.generateDataBase(1.0e-3); + } + else if (sph_system.RestartStep() == 0) + { + write_horizontal_velocity.testResult(); + write_vertical_velocity.testResult(); + } - return 0; + return 0; } diff --git a/tests/2d_examples/test_2d_poiseuille_flow/poiseuille_flow.cpp b/tests/2d_examples/test_2d_poiseuille_flow/poiseuille_flow.cpp index c8b9399169..d4f7a9d2cd 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/poiseuille_flow.cpp +++ b/tests/2d_examples/test_2d_poiseuille_flow/poiseuille_flow.cpp @@ -110,12 +110,14 @@ int main(int ac, char *av[]) Gravity gravity(Vecd(gravity_g, 0.0)); SimpleDynamics constant_gravity(water_block, gravity); SimpleDynamics wall_boundary_normal_direction(wall_boundary); + InteractionWithUpdate kernel_correction_complex(water_block_inner, water_wall_contact); - Dynamics1Level pressure_relaxation(water_block_inner, water_wall_contact); + Dynamics1Level pressure_relaxation(water_block_inner, water_wall_contact); Dynamics1Level density_relaxation(water_block_inner, water_wall_contact); InteractionWithUpdate update_density_by_summation(water_block_inner, water_wall_contact); - InteractionWithUpdate viscous_force(water_block_inner, water_wall_contact); - InteractionWithUpdate> transport_velocity_correction(water_block_inner, water_wall_contact); + InteractionWithUpdate viscous_force(water_block_inner, water_wall_contact); + InteractionWithUpdate> + transport_velocity_correction(water_block_inner, water_wall_contact); ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); ReduceDynamics get_fluid_time_step_size(water_block); @@ -171,6 +173,7 @@ int main(int ac, char *av[]) time_instance = TickCount::now(); Real Dt = get_fluid_advection_time_step_size.exec(); update_density_by_summation.exec(); + kernel_correction_complex.exec(); viscous_force.exec(); transport_velocity_correction.exec(); interval_computing_time_step += TickCount::now() - time_instance; diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml index 780c107a59..21f59ae38f 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml index 5cbc0a88fc..0f3d27c3b9 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml deleted file mode 100644 index b00fe6ed63..0000000000 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml new file mode 100644 index 0000000000..19ffb9658b --- /dev/null +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml index d7aeb54c0d..36b6bf2534 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat index de07de18db..34de7daedb 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat @@ -1,3 +1,3 @@ true -6 +7 4 \ No newline at end of file diff --git a/tests/3d_examples/test_3d_dambreak_elastic_plate_shell/test_3d_dambreak_elastic_plate_shell.cpp b/tests/3d_examples/test_3d_dambreak_elastic_plate_shell/test_3d_dambreak_elastic_plate_shell.cpp index 78a93699c1..10b30403bc 100644 --- a/tests/3d_examples/test_3d_dambreak_elastic_plate_shell/test_3d_dambreak_elastic_plate_shell.cpp +++ b/tests/3d_examples/test_3d_dambreak_elastic_plate_shell/test_3d_dambreak_elastic_plate_shell.cpp @@ -232,7 +232,8 @@ int main(int ac, char *av[]) update_density_by_summation(water_block_inner, water_wall_contact, water_plate_contact); ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); ReduceDynamics get_fluid_time_step_size(water_block); - InteractionWithUpdate, Contact, Contact>, fluid_dynamics::FixedViscosity>> + InteractionWithUpdate, Contact, Contact>, + fluid_dynamics::FixedViscosity, NoKernelCorrection>> viscous_acceleration(water_block_inner, water_wall_contact, water_plate_contact); // FSI InteractionWithUpdate viscous_force_on_plate(plate_water_contact);