From afab322cede30625cf2125ae571963a07c2d2102 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 9 Apr 2024 19:51:37 +0200 Subject: [PATCH 01/10] basic impl --- examples/fsi/falling_rigid_spheres_2d.jl | 118 +++++++++++ src/TrixiParticles.jl | 2 +- src/general/semidiscretization.jl | 99 ++++++--- src/schemes/boundary/rhs.jl | 6 +- src/schemes/schemes.jl | 6 + src/schemes/solid/rigid_body_sph/rhs.jl | 119 +++++++++++ .../solid/rigid_body_sph/rigid_body_sph.jl | 1 + src/schemes/solid/rigid_body_sph/system.jl | 195 ++++++++++++++++++ src/visualization/write2vtk.jl | 8 + 9 files changed, 525 insertions(+), 29 deletions(-) create mode 100644 examples/fsi/falling_rigid_spheres_2d.jl create mode 100644 src/schemes/solid/rigid_body_sph/rhs.jl create mode 100644 src/schemes/solid/rigid_body_sph/rigid_body_sph.jl create mode 100644 src/schemes/solid/rigid_body_sph/system.jl diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl new file mode 100644 index 000000000..a54b06445 --- /dev/null +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -0,0 +1,118 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0075 +solid_particle_spacing = fluid_particle_spacing + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.5) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (2.0, 0.9) +tank_size = (2.0, 1.0) + +fluid_density = 1000.0 +sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere1_radius = 0.3 +sphere2_radius = 0.2 +sphere1_density = 500.0 +sphere2_density = 700.0 + +# Young's modulus and Poisson ratio +sphere1_E = 7e4 +sphere2_E = 1e5 +nu = 0.0 + +sphere1_center = (0.5, 1.6) +sphere2_center = (1.5, 1.6) +sphere1 = SphereShape(solid_particle_spacing, sphere1_radius, sphere1_center, + sphere1_density, sphere_type=VoxelSphere()) +sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, + sphere2_density, sphere_type=VoxelSphere()) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.0 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Solid +solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing +solid_smoothing_kernel = WendlandC2Kernel{2}() + +# For the FSI we need the hydrodynamic masses and densities in the solid boundary model +hydrodynamic_densites_1 = fluid_density * ones(size(sphere1.density)) +hydrodynamic_masses_1 = hydrodynamic_densites_1 * solid_particle_spacing^ndims(fluid_system) + +solid_boundary_model_1 = BoundaryModelDummyParticles(hydrodynamic_densites_1, + hydrodynamic_masses_1, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, + fluid_smoothing_length) + +hydrodynamic_densites_2 = fluid_density * ones(size(sphere2.density)) +hydrodynamic_masses_2 = hydrodynamic_densites_2 * solid_particle_spacing^ndims(fluid_system) + +solid_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2, + hydrodynamic_masses_2, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, + fluid_smoothing_length) + +solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, acceleration=(0.0, -gravity)) +solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL49(), + abstol=1e-6, # Default abstol is 1e-6 + reltol=1e-5, # Default reltol is 1e-3 + save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index b7c0c5fe7..2298cc423 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -42,7 +42,7 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem + BoundarySPHSystem, RigidSPHSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback export ContinuityDensity, SummationDensity diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index ada9275e4..76bfbac60 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -156,7 +156,8 @@ end return compact_support(smoothing_kernel, smoothing_length) end -@inline function compact_support(system::Union{TotalLagrangianSPHSystem, BoundarySPHSystem}, +@inline function compact_support(system::Union{TotalLagrangianSPHSystem, BoundarySPHSystem, + RigidSPHSystem}, neighbor) return compact_support(system, system.boundary_model, neighbor) end @@ -393,6 +394,8 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) + #@trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, du_ode, v_ode, u_ode, semi) + @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) end @@ -470,6 +473,7 @@ function add_source_terms!(dv_ode, v_ode, u_ode, semi) end @inline source_terms(system) = nothing +@inline source_terms(system::RigidSPHSystem) = nothing @inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms @inline add_acceleration!(dv, particle, system) = dv @@ -577,39 +581,82 @@ end end end -# NHS updates -function nhs_coords(system::FluidSystem, - neighbor::FluidSystem, u) - return current_coordinates(u, neighbor) +function collision_interaction!(dv_ode, du_ode, v_ode, u_ode, semi) + # Call `interact!` for each pair of systems + foreach_system(semi) do system + foreach_system(semi) do neighbor + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + if timeit_debug_enabled() + system_index = system_indices(system, semi) + neighbor_index = system_indices(neighbor, semi) + timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" + else + timer_str = "" + end + + collision_interact!(dv_ode, du_ode, v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) + end + end + + return dv_ode end -function nhs_coords(system::FluidSystem, - neighbor::TotalLagrangianSPHSystem, u) - return current_coordinates(u, neighbor) +# FluidSystems don't collide with each other +@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system::FluidSystem, neighbor::FluidSystem, semi; timer_str="") + return dv_ode end -function nhs_coords(system::FluidSystem, - neighbor::BoundarySPHSystem, u) - if neighbor.ismoving[1] - return current_coordinates(u, neighbor) +# BoundarySPHSystems don't collide with each other use RigidSPHSystem +@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system::BoundarySPHSystem, neighbor::BoundarySPHSystem, semi; timer_str="") + return dv_ode +end + + +# Function barrier to make benchmarking interactions easier. +# One can benchmark, e.g. the fluid-fluid interaction, with: +# dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; +# @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); +@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system, neighbor, semi; timer_str="") + dv = wrap_v(dv_ode, system, semi) + du = wrap_u(du_ode, system, semi) + v_system = wrap_v(v_ode, system, semi) + u_system = wrap_u(u_ode, system, semi) + + v_neighbor = wrap_v(v_ode, neighbor, semi) + u_neighbor = wrap_u(u_ode, neighbor, semi) + nhs = get_neighborhood_search(system, neighbor, semi) + + @trixi_timeit timer() timer_str begin + collision_interact!(dv, du, v_system, u_system, v_neighbor, u_neighbor, nhs, system, neighbor) end +end - # Don't update - return nothing +@inline function collision_interact!(dv, du, v_system, u_system, v_neighbor, u_neighbor, nhs, system, neighbor) + return dv end -function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::FluidSystem, u) + +# NHS updates +# All systems that always move update every time +function nhs_coords(system::FluidSystem, + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, RigidSPHSystem}, + u) return current_coordinates(u, neighbor) end function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::TotalLagrangianSPHSystem, u) - # Don't update - return nothing + neighbor::Union{FluidSystem, RigidSPHSystem}, u) + return current_coordinates(u, neighbor) end -function nhs_coords(system::TotalLagrangianSPHSystem, +function nhs_coords(system::RigidSPHSystem, + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, u) + return current_coordinates(u, neighbor) +end + +# Only update when moving +function nhs_coords(system::Union{FluidSystem, TotalLagrangianSPHSystem, RigidSPHSystem}, neighbor::BoundarySPHSystem, u) if neighbor.ismoving[1] return current_coordinates(u, neighbor) @@ -620,7 +667,8 @@ function nhs_coords(system::TotalLagrangianSPHSystem, end function nhs_coords(system::BoundarySPHSystem, - neighbor::FluidSystem, u) + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, + BoundarySPHSystem, RigidSPHSystem}, u) # Don't update return nothing end @@ -630,14 +678,12 @@ function nhs_coords(system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, return current_coordinates(u, neighbor) end -function nhs_coords(system::BoundarySPHSystem, - neighbor::TotalLagrangianSPHSystem, u) +function nhs_coords(system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem, u) # Don't update return nothing end -function nhs_coords(system::BoundarySPHSystem, - neighbor::BoundarySPHSystem, u) +function nhs_coords(system::RigidSPHSystem, neighbor::RigidSPHSystem, u) # Don't update return nothing end @@ -663,7 +709,8 @@ function check_configuration(boundary_system::BoundarySPHSystem, systems) end end -function check_configuration(system::TotalLagrangianSPHSystem, systems) +function check_configuration(system::Union{TotalLagrangianSPHSystem, RigidSPHSystem}, + systems) (; boundary_model) = system foreach_system(systems) do neighbor diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 965e55062..574e832d1 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -1,7 +1,8 @@ # Interaction of boundary with other systems function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::BoundarySPHSystem, neighbor_system) + particle_system::Union{BoundarySPHSystem, RigidSPHSystem}, + neighbor_system) # TODO Solids and moving boundaries should be considered in the continuity equation return dv end @@ -9,7 +10,8 @@ end # For dummy particles with `ContinuityDensity`, solve the continuity equation function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, + particle_system::Union{BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, + RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}}, neighbor_system::FluidSystem) (; boundary_model) = particle_system fluid_density_calculator = neighbor_system.density_calculator diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index 5eabef80e..f7fcc3610 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -2,7 +2,12 @@ # interactions between the different system types. include("fluid/fluid.jl") include("boundary/boundary.jl") + +# needs to be after boundary.jl +include("solid/rigid_body_sph/rigid_body_sph.jl") +# needs to be after boundary.jl include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") + # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem` # and the `TotalLagrangianSPHSystem`. include("boundary/monaghan_kajtar/monaghan_kajtar.jl") @@ -12,3 +17,4 @@ include("fluid/weakly_compressible_sph/rhs.jl") include("fluid/entropically_damped_sph/rhs.jl") include("boundary/rhs.jl") include("solid/total_lagrangian_sph/rhs.jl") +include("solid/rigid_body_sph/rhs.jl") diff --git a/src/schemes/solid/rigid_body_sph/rhs.jl b/src/schemes/solid/rigid_body_sph/rhs.jl new file mode 100644 index 000000000..03a12c675 --- /dev/null +++ b/src/schemes/solid/rigid_body_sph/rhs.jl @@ -0,0 +1,119 @@ + +# Solid-fluid interaction +# Since the object is rigid the movement gets averaged across all particles. +function interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, neighborhood_search, + particle_system::RigidSPHSystem, + neighbor_system::FluidSystem) + sound_speed = system_sound_speed(neighbor_system) + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + total_dv = zeros(ndims(particle_system)) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + for_particle_neighbor(particle_system, neighbor_system, + system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + # Apply the same force to the solid particle + # that the fluid particle experiences due to the soild particle. + # Note that the same arguments are passed here as in fluid-solid interact!, + # except that pos_diff has a flipped sign. + # + # In fluid-solid interaction, use the "hydrodynamic mass" of the solid particles + # corresponding to the rest density of the fluid and not the material density. + m_a = hydrodynamic_mass(particle_system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + rho_a = particle_density(v_particle_system, particle_system, particle) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + + # Use kernel from the fluid system in order to get the same force here in + # solid-fluid interaction as for fluid-solid interaction. + grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance) + + # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles + # corresponding to the chosen boundary model. + p_a = particle_pressure(v_particle_system, particle_system, particle) + p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) + + # Particle and neighbor (and corresponding systems and all corresponding quantities) + # are switched in the following two calls. + # This way, we obtain the exact same force as for the fluid-solid interaction, + # but with a flipped sign (because `pos_diff` is flipped compared to fluid-solid). + dv_boundary = pressure_acceleration(neighbor_system, particle_system, particle, + m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, + distance, grad_kernel, + neighbor_system.correction) + + dv_viscosity_ = dv_viscosity(neighbor_system, particle_system, + v_neighbor_system, v_particle_system, + neighbor, particle, pos_diff, distance, + sound_speed, m_b, m_a, rho_mean) + + dv_particle = dv_boundary + dv_viscosity_ + + for i in 1:ndims(particle_system) + # Multiply `dv` (acceleration on fluid particle b) by the mass of + # particle b to obtain the same force as for the fluid-solid interaction. + # Divide by the material mass of particle a to obtain the acceleration + # of solid particle a. + total_dv[i] += dv[i, particle] + dv_particle[i] * m_b / particle_system.mass[particle] + end + + continuity_equation!(dv, v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system, neighbor_system, grad_kernel) + end + + #total_dv ./=nparticles(particle_system) + for particle in each_moving_particle(particle_system) + for i in 1:ndims(particle_system) + dv[i, particle] += total_dv[i]/nparticles(particle_system) + end + end + return dv +end + +@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::RigidSPHSystem, + neighbor_system::FluidSystem, + grad_kernel) + return dv +end + +@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, + neighbor_system::FluidSystem, + grad_kernel) + fluid_density_calculator = neighbor_system.density_calculator + + v_diff = current_velocity(v_particle_system, particle_system, particle) - + current_velocity(v_neighbor_system, neighbor_system, neighbor) + + # Call the dummy BC version of the continuity equation + continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff, + grad_kernel, particle) +end + +# Solid-boundary interaction +function interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, neighborhood_search, + particle_system::RigidSPHSystem, + neighbor_system::Union{BoundarySPHSystem}) + + + + # TODO continuity equation? + return dv +end diff --git a/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl b/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl new file mode 100644 index 000000000..c67930ab8 --- /dev/null +++ b/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl @@ -0,0 +1 @@ +include("system.jl") diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl new file mode 100644 index 000000000..e4f4b457f --- /dev/null +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -0,0 +1,195 @@ +@doc raw""" + RigidSPHSystem(initial_condition; + boundary_model=nothing, + acceleration=ntuple(_ -> 0.0, NDIMS)) + +System for a SPH particles of a rigid structure. + +# Arguments +- `initial_condition`: Initial condition representing the system's particles. + +# Keyword Arguments +- `boundary_model`: Boundary model to compute the hydrodynamic density and pressure for + fluid-structure interaction (see [Boundary Models](@ref boundary_models)). +- `acceleration`: Acceleration vector for the system. (default: zero vector) +""" +struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} + initial_condition :: InitialCondition{ELTYPE} + local_coordinates :: Array{ELTYPE, 2} # [dimension, particle] + mass :: Array{ELTYPE, 1} # [particle] + material_density :: Array{ELTYPE, 1} # [particle] + acceleration :: SVector{NDIMS, ELTYPE} + center_of_gravity :: SVector{NDIMS, ELTYPE} + boundary_model :: BM + + function RigidSPHSystem(initial_condition; boundary_model=nothing, + acceleration=ntuple(_ -> 0.0, ndims(initial_condition))) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + # Make acceleration an SVector + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end + + local_coordinates = copy(initial_condition.coordinates) + mass = copy(initial_condition.mass) + material_density = copy(initial_condition.density) + + # calculate center of gravity + cog = zeros(SVector{NDIMS, ELTYPE}) + + return new{typeof(boundary_model), NDIMS, ELTYPE}(initial_condition, + local_coordinates, + mass, material_density, + acceleration_, cog, + boundary_model) + end +end + +# Info output functions + +function Base.show(io::IO, system::RigidSPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "RigidSPHSystem{", ndims(system), "}(") + print(io, ", ", system.acceleration) + print(io, ", ", system.boundary_model) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::RigidSPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "RigidSPHSystem{$(ndims(system))}") + summary_line(io, "total #particles", nparticles(system)) + summary_line(io, "acceleration", system.acceleration) + summary_line(io, "boundary model", system.boundary_model) + summary_footer(io) + end +end + +# Memory allocation accessors + +@inline function v_nvariables(system::RigidSPHSystem) + return ndims(system) +end + +@inline function v_nvariables(system::RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}) + return ndims(system) + 1 +end + + +# Variable accessors + +@inline function n_moving_particles(system::RigidSPHSystem) + return nparticles(system) +end + +@inline local_coordinates(system::RigidSPHSystem) = system.local_coordinates + +# @inline function current_coordinates(u, system::RigidSPHSystem) +# return system.coordinates +# end + +@inline function current_velocity(v, system::RigidSPHSystem, particle) + if particle > n_moving_particles(system) + return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) + end + + return extract_svector(v, system, particle) +end + +@inline function viscous_velocity(v, system::RigidSPHSystem, particle) + return extract_svector(system.boundary_model.cache.wall_velocity, system, particle) +end + +@inline function particle_density(v, system::RigidSPHSystem, particle) + return particle_density(v, system.boundary_model, system, particle) +end + +# In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles +# corresponding to the chosen boundary model. +@inline function particle_pressure(v, system::RigidSPHSystem, particle) + return particle_pressure(v, system.boundary_model, system, particle) +end + +@inline function hydrodynamic_mass(system::RigidSPHSystem, particle) + return system.boundary_model.hydrodynamic_mass[particle] +end + +function initialize!(system::RigidSPHSystem, neighborhood_search) + +end + +function write_u0!(u0, system::RigidSPHSystem) + (; initial_condition) = system + + for particle in each_moving_particle(system) + # Write particle coordinates + for dim in 1:ndims(system) + u0[dim, particle] = initial_condition.coordinates[dim, particle] + end + end + + return u0 +end + +function write_v0!(v0, system::RigidSPHSystem) + (; initial_condition, boundary_model) = system + + for particle in each_moving_particle(system) + # Write particle velocities + for dim in 1:ndims(system) + v0[dim, particle] = initial_condition.velocity[dim, particle] + end + end + + write_v0!(v0, boundary_model, system) + + return v0 +end + +function write_v0!(v0, model, system::RigidSPHSystem) + return v0 +end + +function write_v0!(v0, ::BoundaryModelDummyParticles{ContinuityDensity}, + system::RigidSPHSystem) + (; cache) = system.boundary_model + (; initial_density) = cache + + for particle in each_moving_particle(system) + # Set particle densities + v0[ndims(system) + 1, particle] = initial_density[particle] + end + + return v0 +end + +# Update functions +# function update_positions!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) +# # (; movement) = system + +# # movement(system, t) +# end + +# function update_quantities!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) +# return system +# end + +function update_final!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) + (; boundary_model) = system + + # Only update boundary model + update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) +end + + +function viscosity_model(system::RigidSPHSystem) + return system.boundary_model.viscosity +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 9db300269..217544901 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -252,6 +252,14 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end +function write2vtk!(vtk, v, u, t, system::RigidSPHSystem; write_meta_data=true) + + vtk["velocity"] = view(v, 1:ndims(system), :) + vtk["material_density"] = system.material_density + + write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) +end + function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true) write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end From c050291158637a75631dea8aaf1b3e1fd5ec3a6b Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 10 Apr 2024 17:01:43 +0200 Subject: [PATCH 02/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 14 ++- src/general/semidiscretization.jl | 105 ++++++++++++++---- src/schemes/boundary/system.jl | 8 +- .../fluid/weakly_compressible_sph/rhs.jl | 37 ++++++ src/schemes/solid/rigid_body_sph/rhs.jl | 7 +- src/schemes/solid/rigid_body_sph/system.jl | 22 ++-- src/visualization/write2vtk.jl | 1 - 7 files changed, 150 insertions(+), 44 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index a54b06445..523d5f6db 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -3,7 +3,7 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.0075 +fluid_particle_spacing = 0.02 solid_particle_spacing = fluid_particle_spacing # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model @@ -20,7 +20,7 @@ initial_fluid_size = (2.0, 0.9) tank_size = (2.0, 1.0) fluid_density = 1000.0 -sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) +sound_speed = 40 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=1) @@ -32,7 +32,7 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl sphere1_radius = 0.3 sphere2_radius = 0.2 sphere1_density = 500.0 -sphere2_density = 700.0 +sphere2_density = 900.0 # Young's modulus and Poisson ratio sphere1_E = 7e4 @@ -69,7 +69,7 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar boundary_density_calculator, fluid_smoothing_kernel, fluid_smoothing_length) -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Solid @@ -97,8 +97,10 @@ solid_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2, fluid_smoothing_kernel, fluid_smoothing_length) -solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, acceleration=(0.0, -gravity)) -solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, acceleration=(0.0, -gravity)) +solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, + acceleration=(0.0, -gravity), particle_spacing=fluid_particle_spacing) +solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, + acceleration=(0.0, -gravity), particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Simulation diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 76bfbac60..6128136e4 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -394,7 +394,8 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) - #@trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, du_ode, v_ode, u_ode, semi) + @trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, v_ode, + u_ode, semi) @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) end @@ -581,21 +582,26 @@ end end end -function collision_interaction!(dv_ode, du_ode, v_ode, u_ode, semi) +function collision_interaction!(dv_ode, v_ode, u_ode, semi) # Call `interact!` for each pair of systems foreach_system(semi) do system foreach_system(semi) do neighbor - # Construct string for the interactions timer. - # Avoid allocations from string construction when no timers are used. - if timeit_debug_enabled() - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor, semi) - timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" - else - timer_str = "" - end + system_index = system_indices(system, semi) + neighbor_index = system_indices(neighbor, semi) + + # the same system does not collide with its self + if system_index != neighbor_index + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + if timeit_debug_enabled() + timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" + else + timer_str = "" + end - collision_interact!(dv_ode, du_ode, v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) + collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, + timer_str=timer_str) + end end end @@ -603,23 +609,25 @@ function collision_interaction!(dv_ode, du_ode, v_ode, u_ode, semi) end # FluidSystems don't collide with each other -@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system::FluidSystem, neighbor::FluidSystem, semi; timer_str="") +@inline function collision_interact!(dv_ode, v_ode, u_ode, system::FluidSystem, + neighbor::FluidSystem, semi; timer_str="") return dv_ode end # BoundarySPHSystems don't collide with each other use RigidSPHSystem -@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system::BoundarySPHSystem, neighbor::BoundarySPHSystem, semi; timer_str="") +@inline function collision_interact!(dv_ode, v_ode, u_ode, + system::BoundarySPHSystem, neighbor::BoundarySPHSystem, + semi; timer_str="") return dv_ode end - # Function barrier to make benchmarking interactions easier. # One can benchmark, e.g. the fluid-fluid interaction, with: # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; # @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); -@inline function collision_interact!(dv_ode, du_ode, v_ode, u_ode, system, neighbor, semi; timer_str="") +@inline function collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; + timer_str="") dv = wrap_v(dv_ode, system, semi) - du = wrap_u(du_ode, system, semi) v_system = wrap_v(v_ode, system, semi) u_system = wrap_u(u_ode, system, semi) @@ -628,14 +636,73 @@ end nhs = get_neighborhood_search(system, neighbor, semi) @trixi_timeit timer() timer_str begin - collision_interact!(dv, du, v_system, u_system, v_neighbor, u_neighbor, nhs, system, neighbor) + collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, nhs, system, + neighbor) end end -@inline function collision_interact!(dv, du, v_system, u_system, v_neighbor, u_neighbor, nhs, system, neighbor) +@inline function collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, + nhs, system, neighbor) return dv end +# Systems representing solids collide with each other and the boundary +@inline function collision_interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + neighborhood_search, + particle_system::Union{RigidSPHSystem, + TotalLagrangianSPHSystem}, + neighbor_system::Union{BoundarySPHSystem, + RigidSPHSystem, + TotalLagrangianSPHSystem}) + (; particle_spacing) = neighbor_system + neighbor_radius = 0.5 * particle_spacing + + (; particle_spacing) = particle_system + particle_radius = 0.5 * particle_spacing + + collision_distance = neighbor_radius + particle_radius + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + total_dv = zeros(ndims(particle_system)) + + for_particle_neighbor(particle_system, neighbor_system, + system_coords, neighbor_system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + overlap = collision_distance - distance + overlap <= 0 && return + + normal = pos_diff / distance + + part_v = extract_svector(v_particle_system, particle_system, particle) + nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) + rel_vel = part_v - nghbr_v + + # nothing to do we are in contact + norm(rel_vel) <= sqrt(eps()) && return + + println("collide") + + rel_vel_normal = dot(rel_vel, normal) + collision_damping_coefficient = 0.1 + force_magnitude = collision_damping_coefficient * rel_vel_normal + force = force_magnitude * normal + + @inbounds for i in 1:ndims(particle_system) + total_dv[i] += force[i] / mass[particle] + end + end + + for particle in each_moving_particle(particle_system) + for i in 1:ndims(particle_system) + dv[i, particle] += total_dv[i] + end + end + + return dv +end # NHS updates # All systems that always move update every time diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index b30d8fc41..d5a13ef69 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -17,9 +17,10 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM boundary_model :: BM movement :: M ismoving :: Vector{Bool} + particle_spacing :: ELTYPE cache :: C - function BoundarySPHSystem(initial_condition, model; movement=nothing) + function BoundarySPHSystem(initial_condition, model; movement=nothing, particle_spacing = NaN) coordinates = copy(initial_condition.coordinates) NDIMS = size(coordinates, 1) ismoving = zeros(Bool, 1) @@ -33,9 +34,10 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM movement.moving_particles .= collect(1:nparticles(initial_condition)) end + return new{typeof(model), NDIMS, eltype(coordinates), typeof(movement), typeof(cache)}(initial_condition, coordinates, model, movement, - ismoving, cache) + ismoving, particle_spacing, cache) end end @@ -99,6 +101,7 @@ function Base.show(io::IO, system::BoundarySPHSystem) print(io, "BoundarySPHSystem{", ndims(system), "}(") print(io, system.boundary_model) print(io, ", ", system.movement) + print(io, ", ", system.particle_spacing) print(io, ") with ", nparticles(system), " particles") end @@ -111,6 +114,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) summary_header(io, "BoundarySPHSystem{$(ndims(system))}") summary_line(io, "#particles", nparticles(system)) summary_line(io, "boundary model", system.boundary_model) + summary_line(io, "particle spacing", system.particle_spacing) summary_line(io, "movement function", isnothing(system.movement) ? "nothing" : string(system.movement.movement_function)) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0b913bc23..0496d1652 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -153,3 +153,40 @@ end return p_a, p_a end + +# Fluid particles collide with a boundary particle at distance =< 0 to the surface boundary +# particle to prevent them going through the first layer of boundary particles +@inline function collision_interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + neighborhood_search, particle_system::FluidSystem, + neighbor_system::Union{BoundarySPHSystem, + RigidSPHSystem, + TotalLagrangianSPHSystem}) + # (particle_spacing) = neighbor_system + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + for_particle_neighbor(particle_system, neighbor_system, + system_coords, neighbor_system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + if distance <= sqrt(eps()) + # todo: not correct + normal = pos_diff / distance + + part_v = extract_svector(v_particle_system, particle_system, particle) + nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) + # Relative velocity + rel_vel_normal = dot(part_v - nghbr_v, normal) + + collision_damping_coefficient = 0.1 + force_magnitude = collision_damping_coefficient * rel_vel_normal + force = force_magnitude * normal + + @inbounds for i in 1:ndims(particle_system) + dv[i, particle] += force[i] / mass[particle] + end + end + end + return dv +end diff --git a/src/schemes/solid/rigid_body_sph/rhs.jl b/src/schemes/solid/rigid_body_sph/rhs.jl index 03a12c675..acc2bd886 100644 --- a/src/schemes/solid/rigid_body_sph/rhs.jl +++ b/src/schemes/solid/rigid_body_sph/rhs.jl @@ -63,7 +63,8 @@ function interact!(dv, v_particle_system, u_particle_system, # particle b to obtain the same force as for the fluid-solid interaction. # Divide by the material mass of particle a to obtain the acceleration # of solid particle a. - total_dv[i] += dv[i, particle] + dv_particle[i] * m_b / particle_system.mass[particle] + total_dv[i] += dv[i, particle] + + dv_particle[i] * m_b / particle_system.mass[particle] end continuity_equation!(dv, v_particle_system, v_neighbor_system, @@ -75,7 +76,7 @@ function interact!(dv, v_particle_system, u_particle_system, #total_dv ./=nparticles(particle_system) for particle in each_moving_particle(particle_system) for i in 1:ndims(particle_system) - dv[i, particle] += total_dv[i]/nparticles(particle_system) + dv[i, particle] += total_dv[i] / nparticles(particle_system) end end return dv @@ -112,8 +113,6 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::RigidSPHSystem, neighbor_system::Union{BoundarySPHSystem}) - - # TODO continuity equation? return dv end diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl index e4f4b457f..846806952 100644 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -14,16 +14,17 @@ System for a SPH particles of a rigid structure. - `acceleration`: Acceleration vector for the system. (default: zero vector) """ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} - initial_condition :: InitialCondition{ELTYPE} - local_coordinates :: Array{ELTYPE, 2} # [dimension, particle] - mass :: Array{ELTYPE, 1} # [particle] - material_density :: Array{ELTYPE, 1} # [particle] - acceleration :: SVector{NDIMS, ELTYPE} - center_of_gravity :: SVector{NDIMS, ELTYPE} - boundary_model :: BM + initial_condition :: InitialCondition{ELTYPE} + local_coordinates :: Array{ELTYPE, 2} # [dimension, particle] + mass :: Array{ELTYPE, 1} # [particle] + material_density :: Array{ELTYPE, 1} # [particle] + acceleration :: SVector{NDIMS, ELTYPE} + center_of_gravity :: SVector{NDIMS, ELTYPE} + particle_spacing :: ELTYPE + boundary_model :: BM function RigidSPHSystem(initial_condition; boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, ndims(initial_condition))) + acceleration=ntuple(_ -> 0.0, ndims(initial_condition)), particle_spacing=NaN) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -43,7 +44,7 @@ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} return new{typeof(boundary_model), NDIMS, ELTYPE}(initial_condition, local_coordinates, mass, material_density, - acceleration_, cog, + acceleration_, cog, particle_spacing, boundary_model) end end @@ -83,7 +84,6 @@ end return ndims(system) + 1 end - # Variable accessors @inline function n_moving_particles(system::RigidSPHSystem) @@ -123,7 +123,6 @@ end end function initialize!(system::RigidSPHSystem, neighborhood_search) - end function write_u0!(u0, system::RigidSPHSystem) @@ -189,7 +188,6 @@ function update_final!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) end - function viscosity_model(system::RigidSPHSystem) return system.boundary_model.viscosity end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 217544901..b7e7f3f6f 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -253,7 +253,6 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d end function write2vtk!(vtk, v, u, t, system::RigidSPHSystem; write_meta_data=true) - vtk["velocity"] = view(v, 1:ndims(system), :) vtk["material_density"] = system.material_density From 7f2a33ab5d18fbbeaa2ebff0d9f2f130dc854b8e Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 11 Apr 2024 17:06:23 +0200 Subject: [PATCH 03/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 21 ++- src/general/general.jl | 4 + src/general/semidiscretization.jl | 182 +++++++++++++++------ src/schemes/boundary/system.jl | 4 +- src/schemes/solid/rigid_body_sph/system.jl | 10 +- 5 files changed, 161 insertions(+), 60 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index 523d5f6db..fbdef2f92 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -16,7 +16,7 @@ gravity = 9.81 tspan = (0.0, 1.5) # Boundary geometry and initial fluid particle positions -initial_fluid_size = (2.0, 0.9) +initial_fluid_size = (2.0, 0.5) tank_size = (2.0, 1.0) fluid_density = 1000.0 @@ -35,12 +35,12 @@ sphere1_density = 500.0 sphere2_density = 900.0 # Young's modulus and Poisson ratio -sphere1_E = 7e4 -sphere2_E = 1e5 -nu = 0.0 +# sphere1_E = 7e4 +# sphere2_E = 1e5 +# nu = 0.0 -sphere1_center = (0.5, 1.6) -sphere2_center = (1.5, 1.6) +sphere1_center = (0.5, 1.1) +sphere2_center = (1.5, 1.1) sphere1 = SphereShape(solid_particle_spacing, sphere1_radius, sphere1_center, sphere1_density, sphere_type=VoxelSphere()) sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, @@ -69,7 +69,8 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar boundary_density_calculator, fluid_smoothing_kernel, fluid_smoothing_length) -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, particle_spacing=fluid_particle_spacing) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Solid @@ -98,9 +99,11 @@ solid_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2, fluid_smoothing_length) solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, - acceleration=(0.0, -gravity), particle_spacing=fluid_particle_spacing) + acceleration=(0.0, -gravity), + particle_spacing=fluid_particle_spacing) solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, - acceleration=(0.0, -gravity), particle_spacing=fluid_particle_spacing) + acceleration=(0.0, -gravity), + particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Simulation diff --git a/src/general/general.jl b/src/general/general.jl index 6651c03c5..9576cf2c7 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -15,6 +15,10 @@ timer_name(::BoundarySystem) = "boundary" return du end +mutable struct MutableBool + value::Bool +end + # Note that `semidiscretization.jl` depends on the system types and has to be # included later. # `density_calculators.jl` needs to be included before `corrections.jl`. diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6128136e4..e68954ad7 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -49,6 +49,7 @@ struct Semidiscretization{S, RU, RV, NS} ranges_u :: RU ranges_v :: RV neighborhood_searches :: NS + systems_have_collided :: MutableBool # Dispatch at `systems` to distinguish this constructor from the one below when # 4 systems are passed. @@ -56,7 +57,7 @@ struct Semidiscretization{S, RU, RV, NS} function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches) new{typeof(systems), typeof(ranges_u), typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v, - neighborhood_searches) + neighborhood_searches, MutableBool(false)) end end @@ -394,10 +395,10 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) - @trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, v_ode, - u_ode, semi) - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) + + @trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, v_ode, + u_ode, semi) end return dv_ode @@ -582,32 +583,58 @@ end end end +# start of recursion function collision_interaction!(dv_ode, v_ode, u_ode, semi) - # Call `interact!` for each pair of systems + + collision_interaction!(dv_ode, v_ode, u_ode, semi, true) + + return dv_ode +end + +@inline function reset_system_collision_status(semi) foreach_system(semi) do system - foreach_system(semi) do neighbor - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor, semi) - - # the same system does not collide with its self - if system_index != neighbor_index - # Construct string for the interactions timer. - # Avoid allocations from string construction when no timers are used. - if timeit_debug_enabled() - timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" - else - timer_str = "" - end + if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem + system.has_collided.value = false + end + end +end + +function collision_interaction!(dv_ode, v_ode, u_ode, semi, systems_have_collided) + if systems_have_collided + semi.systems_have_collided.value = false + reset_system_collision_status(semi) - collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, - timer_str=timer_str) + # Call `interact!` for each pair of systems + foreach_system(semi) do system + foreach_system(semi) do neighbor + system_index = system_indices(system, semi) + neighbor_index = system_indices(neighbor, semi) + + # the same system does not collide with its self + if system_index != neighbor_index + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + if timeit_debug_enabled() + timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" + else + timer_str = "" + end + + collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, + timer_str=timer_str) + end end end + + # if systems have collided we need to call this again since there can be more collisions + collision_interaction!(dv_ode, v_ode, u_ode, semi, semi.systems_have_collided.value) end return dv_ode end + + # FluidSystems don't collide with each other @inline function collision_interact!(dv_ode, v_ode, u_ode, system::FluidSystem, neighbor::FluidSystem, semi; timer_str="") @@ -621,6 +648,15 @@ end return dv_ode end +@inline function set_collision_status!(semi, system) + return semi +end + +@inline function set_collision_status!(semi, system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}) + semi.systems_have_collided.value = semi.systems_have_collided.value || system.has_collided.value + return semi +end + # Function barrier to make benchmarking interactions easier. # One can benchmark, e.g. the fluid-fluid interaction, with: # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; @@ -639,6 +675,8 @@ end collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, nhs, system, neighbor) end + + set_collision_status!(semi, system) end @inline function collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, @@ -646,15 +684,13 @@ end return dv end -# Systems representing solids collide with each other and the boundary +# Systems representing solids colliding with a boundary under the assumption the boundary is not influenced by the collision @inline function collision_interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}, - neighbor_system::Union{BoundarySPHSystem, - RigidSPHSystem, - TotalLagrangianSPHSystem}) + neighbor_system::BoundarySPHSystem) (; particle_spacing) = neighbor_system neighbor_radius = 0.5 * particle_spacing @@ -668,38 +704,92 @@ end total_dv = zeros(ndims(particle_system)) - for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_system_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance + max_overlap = 0.0 + max_normal = zeros(ndims(particle_system)) + + # Detect the maximum overlap and calculate the normal vector at maximum overlap + for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance overlap = collision_distance - distance - overlap <= 0 && return + if overlap > max_overlap + max_overlap = overlap + max_normal = pos_diff / distance # Normal vector pointing from particle to neighbor + end + end - normal = pos_diff / distance + # Only proceed if there is an overlap + if max_overlap > 0 + # current_dv_normal = dot(dv[:, particle], max_normal) * max_normal - part_v = extract_svector(v_particle_system, particle_system, particle) - nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) - rel_vel = part_v - nghbr_v + # Calculate the total initial velocity in the direction of the normal + particle_velocity = extract_svector(v_particle_system, particle_system, 1) + initial_normal_velocity = dot(particle_velocity, max_normal) - # nothing to do we are in contact - norm(rel_vel) <= sqrt(eps()) && return + # Coefficient of restitution (0 < e <= 1) + # e = 1 -> perfectly elastic collision + # 0 < e < 1 -> inelastic collision - println("collide") + e = 1.0 + if norm(particle_velocity) < sqrt(eps()) + e= 0.0 + end - rel_vel_normal = dot(rel_vel, normal) - collision_damping_coefficient = 0.1 - force_magnitude = collision_damping_coefficient * rel_vel_normal - force = force_magnitude * normal + # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) + change_in_velocity_normal = -1000 * (1 + e) * initial_normal_velocity * max_normal - @inbounds for i in 1:ndims(particle_system) - total_dv[i] += force[i] / mass[particle] + # Apply the change uniformly across all particles + @inbounds for particle in each_moving_particle(particle_system) + dv[:, particle] = change_in_velocity_normal end end - for particle in each_moving_particle(particle_system) - for i in 1:ndims(particle_system) - dv[i, particle] += total_dv[i] - end - end + + + # doesn't work + # for particle in each_moving_particle(particle_system) + # @inbounds for i in 1:ndims(particle_system) + # u_particle_system[i, particle] -= position_correction[i] + # end + # end + + + + # normal = pos_diff / distance + + # part_v = extract_svector(v_particle_system, particle_system, particle) + # nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) + # rel_vel = part_v - nghbr_v + + # # nothing to do we are in contact + # norm(rel_vel) <= sqrt(eps()) && return + + # println("collide") + + # rel_vel_normal = dot(rel_vel, normal) + # collision_damping_coefficient = 0.1 + # force_magnitude = collision_damping_coefficient * rel_vel_normal + # force = force_magnitude * normal + + # @inbounds for i in 1:ndims(particle_system) + # total_dv[i] += force[i] / mass[particle] + # end + + # for particle in each_moving_particle(particle_system) + # for i in 1:ndims(particle_system) + # dv[i, particle] += total_dv[i] + # end + # end + + return dv +end + +# collision influences both systems +@inline function collision_interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + neighborhood_search, + particle_system::Union{RigidSPHSystem, + TotalLagrangianSPHSystem}, + neighbor_system::Union{RigidSPHSystem, + TotalLagrangianSPHSystem}) return dv end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index d5a13ef69..1e5a325aa 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -20,7 +20,8 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM particle_spacing :: ELTYPE cache :: C - function BoundarySPHSystem(initial_condition, model; movement=nothing, particle_spacing = NaN) + function BoundarySPHSystem(initial_condition, model; movement=nothing, + particle_spacing=NaN) coordinates = copy(initial_condition.coordinates) NDIMS = size(coordinates, 1) ismoving = zeros(Bool, 1) @@ -34,7 +35,6 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM movement.moving_particles .= collect(1:nparticles(initial_condition)) end - return new{typeof(model), NDIMS, eltype(coordinates), typeof(movement), typeof(cache)}(initial_condition, coordinates, model, movement, ismoving, particle_spacing, cache) diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl index 846806952..3baa6e137 100644 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -21,10 +21,12 @@ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} acceleration :: SVector{NDIMS, ELTYPE} center_of_gravity :: SVector{NDIMS, ELTYPE} particle_spacing :: ELTYPE + has_collided :: MutableBool boundary_model :: BM function RigidSPHSystem(initial_condition; boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, ndims(initial_condition)), particle_spacing=NaN) + acceleration=ntuple(_ -> 0.0, ndims(initial_condition)), + particle_spacing=NaN) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -38,13 +40,15 @@ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} mass = copy(initial_condition.mass) material_density = copy(initial_condition.density) - # calculate center of gravity + # TODO: calculate center of gravity cog = zeros(SVector{NDIMS, ELTYPE}) return new{typeof(boundary_model), NDIMS, ELTYPE}(initial_condition, local_coordinates, mass, material_density, - acceleration_, cog, particle_spacing, + acceleration_, cog, + particle_spacing, + MutableBool(false), boundary_model) end end From 6224f6c668e2df0eb8723e572181195cb485e9a4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 12 Apr 2024 08:11:47 +0200 Subject: [PATCH 04/10] update --- src/general/semidiscretization.jl | 11 ++- .../fluid/weakly_compressible_sph/rhs.jl | 72 +++++++++---------- 2 files changed, 46 insertions(+), 37 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index e68954ad7..e0a51b4ee 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -370,6 +370,8 @@ function drift!(du_ode, v_ode, u_ode, semi, t) end end end + # @trixi_timeit timer() "collision interaction" collision_interaction!(du_ode, v_ode, + # u_ode, semi) end return du_ode @@ -664,6 +666,8 @@ end @inline function collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; timer_str="") dv = wrap_v(dv_ode, system, semi) + # dv = wrap_u(dv_ode, system, semi) + v_system = wrap_v(v_ode, system, semi) u_system = wrap_u(u_ode, system, semi) @@ -734,10 +738,15 @@ end end # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) - change_in_velocity_normal = -1000 * (1 + e) * initial_normal_velocity * max_normal + change_in_velocity_normal = -(1 + e) * initial_normal_velocity * max_normal + + #collision_time = max_overlap/norm(particle_velocity) + #println("collison_time", collision_time) # Apply the change uniformly across all particles @inbounds for particle in each_moving_particle(particle_system) + # current_dv_normal = dot(dv[:, particle], max_normal) * max_normal + # dv[:, particle] = change_in_velocity_normal/collision_time dv[:, particle] = change_in_velocity_normal end end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0496d1652..21ed2d747 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -154,39 +154,39 @@ end return p_a, p_a end -# Fluid particles collide with a boundary particle at distance =< 0 to the surface boundary -# particle to prevent them going through the first layer of boundary particles -@inline function collision_interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search, particle_system::FluidSystem, - neighbor_system::Union{BoundarySPHSystem, - RigidSPHSystem, - TotalLagrangianSPHSystem}) - # (particle_spacing) = neighbor_system - - system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - - for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_system_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - if distance <= sqrt(eps()) - # todo: not correct - normal = pos_diff / distance - - part_v = extract_svector(v_particle_system, particle_system, particle) - nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) - # Relative velocity - rel_vel_normal = dot(part_v - nghbr_v, normal) - - collision_damping_coefficient = 0.1 - force_magnitude = collision_damping_coefficient * rel_vel_normal - force = force_magnitude * normal - - @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += force[i] / mass[particle] - end - end - end - return dv -end +# # Fluid particles collide with a boundary particle at distance =< 0 to the surface boundary +# # particle to prevent them going through the first layer of boundary particles +# @inline function collision_interact!(dv, v_particle_system, u_particle_system, +# v_neighbor_system, u_neighbor_system, +# neighborhood_search, particle_system::FluidSystem, +# neighbor_system::Union{BoundarySPHSystem, +# RigidSPHSystem, +# TotalLagrangianSPHSystem}) +# # (particle_spacing) = neighbor_system + +# system_coords = current_coordinates(u_particle_system, particle_system) +# neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + +# for_particle_neighbor(particle_system, neighbor_system, +# system_coords, neighbor_system_coords, +# neighborhood_search) do particle, neighbor, pos_diff, distance +# if distance <= sqrt(eps()) +# # todo: not correct +# normal = pos_diff / distance + +# part_v = extract_svector(v_particle_system, particle_system, particle) +# nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) +# # Relative velocity +# rel_vel_normal = dot(part_v - nghbr_v, normal) + +# collision_damping_coefficient = 0.1 +# force_magnitude = collision_damping_coefficient * rel_vel_normal +# force = force_magnitude * normal + +# @inbounds for i in 1:ndims(particle_system) +# dv[i, particle] += force[i] / mass[particle] +# end +# end +# end +# return dv +# end From a92235f17689aa0662f399da3390480a2ddf6291 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 12 Apr 2024 16:31:40 +0200 Subject: [PATCH 05/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 52 ++++++++----- src/general/semidiscretization.jl | 86 +++++++++++++--------- src/schemes/solid/rigid_body_sph/system.jl | 8 +- 3 files changed, 90 insertions(+), 56 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index fbdef2f92..c6bc27639 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -3,7 +3,7 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.02 +fluid_particle_spacing = 0.0125 solid_particle_spacing = fluid_particle_spacing # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model @@ -13,14 +13,14 @@ spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -tspan = (0.0, 1.5) +tspan = (0.0, 2.0) # Boundary geometry and initial fluid particle positions initial_fluid_size = (2.0, 0.5) -tank_size = (2.0, 1.0) +tank_size = (2.0, 2.0) fluid_density = 1000.0 -sound_speed = 40 +sound_speed = 100 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=1) @@ -31,13 +31,11 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl sphere1_radius = 0.3 sphere2_radius = 0.2 -sphere1_density = 500.0 -sphere2_density = 900.0 - -# Young's modulus and Poisson ratio -# sphere1_E = 7e4 -# sphere2_E = 1e5 -# nu = 0.0 +# wood +sphere1_density = 600.0 +# steel +#sphere2_density = 7700 +sphere2_density = 1500 sphere1_center = (0.5, 1.1) sphere2_center = (1.5, 1.1) @@ -52,7 +50,7 @@ fluid_smoothing_length = 3.0 * fluid_particle_spacing fluid_smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +viscosity = ViscosityAdami(nu=1e-4) density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, @@ -67,15 +65,13 @@ boundary_density_calculator = AdamiPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length) + fluid_smoothing_kernel, fluid_smoothing_length, viscosity=viscosity) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Solid -solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing -solid_smoothing_kernel = WendlandC2Kernel{2}() # For the FSI we need the hydrodynamic masses and densities in the solid boundary model hydrodynamic_densites_1 = fluid_density * ones(size(sphere1.density)) @@ -86,7 +82,7 @@ solid_boundary_model_1 = BoundaryModelDummyParticles(hydrodynamic_densites_1, state_equation=state_equation, boundary_density_calculator, fluid_smoothing_kernel, - fluid_smoothing_length) + fluid_smoothing_length, viscosity=viscosity) hydrodynamic_densites_2 = fluid_density * ones(size(sphere2.density)) hydrodynamic_masses_2 = hydrodynamic_densites_2 * solid_particle_spacing^ndims(fluid_system) @@ -96,7 +92,7 @@ solid_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2, state_equation=state_equation, boundary_density_calculator, fluid_smoothing_kernel, - fluid_smoothing_length) + fluid_smoothing_length,viscosity=viscosity) solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, acceleration=(0.0, -gravity), @@ -116,8 +112,26 @@ saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix callbacks = CallbackSet(info_callback, saving_callback) +function collision(vu, integrator, semi, t) + v_ode, u_ode = vu.x + + TrixiParticles.foreach_system(semi) do system + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + + if system isa RigidSPHSystem && system.has_collided.value + println(system.collision_u) + + for particle in TrixiParticles.each_moving_particle(system) + v[:, particle] += system.collision_impulse + u[:, particle] += system.collision_u + end + end + end +end + # Use a Runge-Kutta method with automatic (error based) time step size control. -sol = solve(ode, RDPK3SpFSAL49(), +sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), abstol=1e-6, # Default abstol is 1e-6 - reltol=1e-5, # Default reltol is 1e-3 + reltol=1e-4, # Default reltol is 1e-3 save_everystep=false, callback=callbacks); diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index e0a51b4ee..709acd258 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -57,7 +57,8 @@ struct Semidiscretization{S, RU, RV, NS} function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches) new{typeof(systems), typeof(ranges_u), typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v, - neighborhood_searches, MutableBool(false)) + neighborhood_searches, + MutableBool(false)) end end @@ -354,6 +355,10 @@ function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) end +function calculate_dt(v_ode, u_ode, cfl_number, system) + return Inf +end + function drift!(du_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "drift!" begin @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode) @@ -368,8 +373,14 @@ function drift!(du_ode, v_ode, u_ode, semi, t) # This can be dispatched per system add_velocity!(du, v, particle, system) end + # if system isa RigidSPHSystem && system.has_collided.value + # for particle in each_moving_particle(system) + # du[:, particle] = system.collision_impulse + # end + # end end end + # @trixi_timeit timer() "collision interaction" collision_interaction!(du_ode, v_ode, # u_ode, semi) end @@ -400,7 +411,7 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) @trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, v_ode, - u_ode, semi) + u_ode, semi) end return dv_ode @@ -587,7 +598,6 @@ end # start of recursion function collision_interaction!(dv_ode, v_ode, u_ode, semi) - collision_interaction!(dv_ode, v_ode, u_ode, semi, true) return dv_ode @@ -635,8 +645,6 @@ function collision_interaction!(dv_ode, v_ode, u_ode, semi, systems_have_collide return dv_ode end - - # FluidSystems don't collide with each other @inline function collision_interact!(dv_ode, v_ode, u_ode, system::FluidSystem, neighbor::FluidSystem, semi; timer_str="") @@ -650,14 +658,14 @@ end return dv_ode end -@inline function set_collision_status!(semi, system) - return semi -end +# @inline function set_collision_status!(semi, system) +# return semi +# end -@inline function set_collision_status!(semi, system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}) - semi.systems_have_collided.value = semi.systems_have_collided.value || system.has_collided.value - return semi -end +# @inline function set_collision_status!(semi, system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}) +# semi.systems_have_collided.value = semi.systems_have_collided.value || system.has_collided.value +# return semi +# end # Function barrier to make benchmarking interactions easier. # One can benchmark, e.g. the fluid-fluid interaction, with: @@ -680,7 +688,7 @@ end neighbor) end - set_collision_status!(semi, system) + # set_collision_status!(semi, system) end @inline function collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, @@ -710,9 +718,13 @@ end max_overlap = 0.0 max_normal = zeros(ndims(particle_system)) + particle_system.collision_impulse .= zeros(ndims(particle_system)) + particle_system.collision_u .= zeros(ndims(particle_system)) # Detect the maximum overlap and calculate the normal vector at maximum overlap - for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance + for_particle_neighbor(particle_system, neighbor_system, system_coords, + neighbor_system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance overlap = collision_distance - distance if overlap > max_overlap max_overlap = overlap @@ -720,21 +732,25 @@ end end end - # Only proceed if there is an overlap if max_overlap > 0 + particle_system.has_collided.value = true # current_dv_normal = dot(dv[:, particle], max_normal) * max_normal - # Calculate the total initial velocity in the direction of the normal particle_velocity = extract_svector(v_particle_system, particle_system, 1) initial_normal_velocity = dot(particle_velocity, max_normal) - # Coefficient of restitution (0 < e <= 1) - # e = 1 -> perfectly elastic collision - # 0 < e < 1 -> inelastic collision + # we need to move in the opposite direction the same amount so times 2 + particle_system.collision_u .= (max_overlap + sqrt(eps()))* max_normal + max_overlap/norm(particle_velocity) * initial_normal_velocity - e = 1.0 + # Coefficient of restitution (0 < e <= 1) + # Tungsten Carbide = 0.7 to 1.0 + # Steel = 0.6 to 0.9 + # High Grade aluminium hardened = 0.7 to 0.8 + # Low Grade aluminium unhardened = 0.3 to 0.4 + # Pure aluminium = 0.08 to 0.12 + e = 0.8 if norm(particle_velocity) < sqrt(eps()) - e= 0.0 + e = 0.0 end # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) @@ -743,16 +759,17 @@ end #collision_time = max_overlap/norm(particle_velocity) #println("collison_time", collision_time) + particle_system.collision_impulse .= change_in_velocity_normal + # Apply the change uniformly across all particles - @inbounds for particle in each_moving_particle(particle_system) - # current_dv_normal = dot(dv[:, particle], max_normal) * max_normal - # dv[:, particle] = change_in_velocity_normal/collision_time - dv[:, particle] = change_in_velocity_normal + for particle in each_moving_particle(particle_system) + current_dv_normal = dot(dv[:, particle], max_normal) * max_normal + # dv[:, particle] = change_in_velocity_normal / collision_time + # dv[:, particle] += -current_dv_normal + change_in_velocity_normal * 50000 + dv[:, particle] += -current_dv_normal end end - - # doesn't work # for particle in each_moving_particle(particle_system) # @inbounds for i in 1:ndims(particle_system) @@ -760,8 +777,6 @@ end # end # end - - # normal = pos_diff / distance # part_v = extract_svector(v_particle_system, particle_system, particle) @@ -793,13 +808,12 @@ end # collision influences both systems @inline function collision_interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search, - particle_system::Union{RigidSPHSystem, - TotalLagrangianSPHSystem}, - neighbor_system::Union{RigidSPHSystem, - TotalLagrangianSPHSystem}) - + v_neighbor_system, u_neighbor_system, + neighborhood_search, + particle_system::Union{RigidSPHSystem, + TotalLagrangianSPHSystem}, + neighbor_system::Union{RigidSPHSystem, + TotalLagrangianSPHSystem}) return dv end diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl index 3baa6e137..9ff0451fd 100644 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -19,7 +19,9 @@ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} mass :: Array{ELTYPE, 1} # [particle] material_density :: Array{ELTYPE, 1} # [particle] acceleration :: SVector{NDIMS, ELTYPE} - center_of_gravity :: SVector{NDIMS, ELTYPE} + center_of_gravity :: Array{ELTYPE, 1} # [dimension] + collision_impulse :: Array{ELTYPE, 1} # [dimension] + collision_u :: Array{ELTYPE, 1} # [dimension] particle_spacing :: ELTYPE has_collided :: MutableBool boundary_model :: BM @@ -42,11 +44,15 @@ struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} # TODO: calculate center of gravity cog = zeros(SVector{NDIMS, ELTYPE}) + collision_impulse = zeros(SVector{NDIMS, ELTYPE}) + collision_u = zeros(SVector{NDIMS, ELTYPE}) return new{typeof(boundary_model), NDIMS, ELTYPE}(initial_condition, local_coordinates, mass, material_density, acceleration_, cog, + collision_impulse, + collision_u, particle_spacing, MutableBool(false), boundary_model) From 32969d0db9a5149c07b05649b14bef77ef505f5a Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 15 Apr 2024 13:16:53 +0200 Subject: [PATCH 06/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 18 +++--- src/general/semidiscretization.jl | 73 ++++++------------------ 2 files changed, 26 insertions(+), 65 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index c6bc27639..4dfd7a749 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -3,7 +3,7 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.0125 +fluid_particle_spacing = 0.015 solid_particle_spacing = fluid_particle_spacing # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model @@ -13,14 +13,14 @@ spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -tspan = (0.0, 2.0) +tspan = (0.0, 4.0) # Boundary geometry and initial fluid particle positions initial_fluid_size = (2.0, 0.5) tank_size = (2.0, 2.0) fluid_density = 1000.0 -sound_speed = 100 +sound_speed = 150 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=1) @@ -35,10 +35,10 @@ sphere2_radius = 0.2 sphere1_density = 600.0 # steel #sphere2_density = 7700 -sphere2_density = 1500 +sphere2_density = 3000 sphere1_center = (0.5, 1.1) -sphere2_center = (1.5, 1.1) +sphere2_center = (1.5, 0.8) sphere1 = SphereShape(solid_particle_spacing, sphere1_radius, sphere1_center, sphere1_density, sphere_type=VoxelSphere()) sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, @@ -103,7 +103,7 @@ solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2) +semi = Semidiscretization(boundary_system, solid_system_2) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) @@ -120,8 +120,6 @@ function collision(vu, integrator, semi, t) u = TrixiParticles.wrap_u(u_ode, system, semi) if system isa RigidSPHSystem && system.has_collided.value - println(system.collision_u) - for particle in TrixiParticles.each_moving_particle(system) v[:, particle] += system.collision_impulse u[:, particle] += system.collision_u @@ -132,6 +130,6 @@ end # Use a Runge-Kutta method with automatic (error based) time step size control. sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), - abstol=1e-6, # Default abstol is 1e-6 - reltol=1e-4, # Default reltol is 1e-3 + abstol=1e-7, # Default abstol is 1e-6 + reltol=1e-5, # Default reltol is 1e-3 save_everystep=false, callback=callbacks); diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 709acd258..bab9be5bc 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -714,12 +714,10 @@ end system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - total_dv = zeros(ndims(particle_system)) - - max_overlap = 0.0 - max_normal = zeros(ndims(particle_system)) particle_system.collision_impulse .= zeros(ndims(particle_system)) particle_system.collision_u .= zeros(ndims(particle_system)) + max_overlap = 0.0 + collision_normal = zeros(ndims(particle_system)) # Detect the maximum overlap and calculate the normal vector at maximum overlap for_particle_neighbor(particle_system, neighbor_system, system_coords, @@ -728,19 +726,15 @@ end overlap = collision_distance - distance if overlap > max_overlap max_overlap = overlap - max_normal = pos_diff / distance # Normal vector pointing from particle to neighbor + # Normal vector pointing from particle to neighbor + collision_normal = pos_diff / distance end end if max_overlap > 0 particle_system.has_collided.value = true - # current_dv_normal = dot(dv[:, particle], max_normal) * max_normal - particle_velocity = extract_svector(v_particle_system, particle_system, 1) - initial_normal_velocity = dot(particle_velocity, max_normal) - - # we need to move in the opposite direction the same amount so times 2 - particle_system.collision_u .= (max_overlap + sqrt(eps()))* max_normal + max_overlap/norm(particle_velocity) * initial_normal_velocity + initial_velocity_normal = dot(extract_svector(v_particle_system, particle_system, 1), collision_normal) # Coefficient of restitution (0 < e <= 1) # Tungsten Carbide = 0.7 to 1.0 @@ -748,61 +742,30 @@ end # High Grade aluminium hardened = 0.7 to 0.8 # Low Grade aluminium unhardened = 0.3 to 0.4 # Pure aluminium = 0.08 to 0.12 - e = 0.8 - if norm(particle_velocity) < sqrt(eps()) - e = 0.0 + restitution_coefficient = 0.1 + if norm(initial_velocity_normal) < 10*sqrt(eps()) + restitution_coefficient = 0.0 end # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) - change_in_velocity_normal = -(1 + e) * initial_normal_velocity * max_normal + velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * collision_normal - #collision_time = max_overlap/norm(particle_velocity) - #println("collison_time", collision_time) + # we need to move in the opposite direction the same amount so times 2 + positional_correction = (2.0 * max_overlap + sqrt(eps())) * collision_normal + # if we don't move back far enough the current version gets unstable + movement_scale = max(1.0, 0.25 * particle_radius / norm(positional_correction)) + positional_correction_scaled = positional_correction * movement_scale - particle_system.collision_impulse .= change_in_velocity_normal + particle_system.collision_u .= positional_correction_scaled + particle_system.collision_impulse .= velocity_change_normal # Apply the change uniformly across all particles for particle in each_moving_particle(particle_system) - current_dv_normal = dot(dv[:, particle], max_normal) * max_normal - # dv[:, particle] = change_in_velocity_normal / collision_time - # dv[:, particle] += -current_dv_normal + change_in_velocity_normal * 50000 - dv[:, particle] += -current_dv_normal + current_normal_component = dot(dv[:, particle], collision_normal) * collision_normal + dv[:, particle] += -current_normal_component end end - # doesn't work - # for particle in each_moving_particle(particle_system) - # @inbounds for i in 1:ndims(particle_system) - # u_particle_system[i, particle] -= position_correction[i] - # end - # end - - # normal = pos_diff / distance - - # part_v = extract_svector(v_particle_system, particle_system, particle) - # nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) - # rel_vel = part_v - nghbr_v - - # # nothing to do we are in contact - # norm(rel_vel) <= sqrt(eps()) && return - - # println("collide") - - # rel_vel_normal = dot(rel_vel, normal) - # collision_damping_coefficient = 0.1 - # force_magnitude = collision_damping_coefficient * rel_vel_normal - # force = force_magnitude * normal - - # @inbounds for i in 1:ndims(particle_system) - # total_dv[i] += force[i] / mass[particle] - # end - - # for particle in each_moving_particle(particle_system) - # for i in 1:ndims(particle_system) - # dv[i, particle] += total_dv[i] - # end - # end - return dv end From adc9f5e8103e1f3461ac97d195682fefb4a723e0 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 17 Apr 2024 01:56:11 +0200 Subject: [PATCH 07/10] debugging --- examples/fsi/falling_rigid_spheres_2d.jl | 39 ++-- src/general/semidiscretization.jl | 214 +++++++++++++++------ src/schemes/solid/rigid_body_sph/system.jl | 4 - 3 files changed, 181 insertions(+), 76 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index 4dfd7a749..2c68ce3e4 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -1,5 +1,6 @@ using TrixiParticles using OrdinaryDiffEq +using LinearAlgebra # ========================================================================================== # ==== Resolution @@ -13,7 +14,7 @@ spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -tspan = (0.0, 4.0) +tspan = (0.0, 0.4) # Boundary geometry and initial fluid particle positions initial_fluid_size = (2.0, 0.5) @@ -107,7 +108,7 @@ semi = Semidiscretization(boundary_system, solid_system_2) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) -saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="", +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", write_meta_data=true) callbacks = CallbackSet(info_callback, saving_callback) @@ -115,17 +116,29 @@ callbacks = CallbackSet(info_callback, saving_callback) function collision(vu, integrator, semi, t) v_ode, u_ode = vu.x - TrixiParticles.foreach_system(semi) do system - v = TrixiParticles.wrap_v(v_ode, system, semi) - u = TrixiParticles.wrap_u(u_ode, system, semi) - - if system isa RigidSPHSystem && system.has_collided.value - for particle in TrixiParticles.each_moving_particle(system) - v[:, particle] += system.collision_impulse - u[:, particle] += system.collision_u - end - end - end + println("vollision at ", t) + + TrixiParticles.collision_interaction!(v_ode, u_ode, semi) + + # TrixiParticles.foreach_system(semi) do system + # v = TrixiParticles.wrap_v(v_ode, system, semi) + # u = TrixiParticles.wrap_u(u_ode, system, semi) + + # if system isa RigidSPHSystem && system.has_collided.value + # velocity_change = norm(v[:, 1]) - norm(v[:, 1] + system.collision_impulse) + # if abs(velocity_change) > 1 + # println("before: ", v[:, 1]) + # println("after: ", v[:, 1] + system.collision_impulse) + # println("imp: ", system.collision_impulse) + + # exit(-1) + # end + # for particle in TrixiParticles.each_moving_particle(system) + # v[:, particle] += system.collision_impulse + # u[:, particle] += system.collision_u + # end + # end + # end end # Use a Runge-Kutta method with automatic (error based) time step size control. diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index bab9be5bc..b236c7b95 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -373,11 +373,9 @@ function drift!(du_ode, v_ode, u_ode, semi, t) # This can be dispatched per system add_velocity!(du, v, particle, system) end - # if system isa RigidSPHSystem && system.has_collided.value - # for particle in each_moving_particle(system) - # du[:, particle] = system.collision_impulse - # end - # end + if system isa RigidSPHSystem && system.has_collided.value + println("body speed", v[:,1]) + end end end @@ -399,6 +397,7 @@ end @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du function kick!(dv_ode, v_ode, u_ode, semi, t) + println("kick at ", t) @trixi_timeit timer() "kick!" begin @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode) @@ -409,9 +408,6 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) semi) @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) - - @trixi_timeit timer() "collision interaction" collision_interaction!(dv_ode, v_ode, - u_ode, semi) end return dv_ode @@ -597,24 +593,24 @@ end end # start of recursion -function collision_interaction!(dv_ode, v_ode, u_ode, semi) - collision_interaction!(dv_ode, v_ode, u_ode, semi, true) +function collision_interaction!(v_ode, u_ode, semi) + collision_interaction!(v_ode, u_ode, semi, true) - return dv_ode + return v_ode end -@inline function reset_system_collision_status(semi) - foreach_system(semi) do system - if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem - system.has_collided.value = false - end - end -end +# @inline function reset_system_collision_status(semi) +# foreach_system(semi) do system +# if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem +# system.has_collided.value = false +# end +# end +# end -function collision_interaction!(dv_ode, v_ode, u_ode, semi, systems_have_collided) +function collision_interaction!(v_ode, u_ode, semi, systems_have_collided) if systems_have_collided semi.systems_have_collided.value = false - reset_system_collision_status(semi) + # reset_system_collision_status(semi) # Call `interact!` for each pair of systems foreach_system(semi) do system @@ -632,30 +628,30 @@ function collision_interaction!(dv_ode, v_ode, u_ode, semi, systems_have_collide timer_str = "" end - collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, + collision_interact!(v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) end end end # if systems have collided we need to call this again since there can be more collisions - collision_interaction!(dv_ode, v_ode, u_ode, semi, semi.systems_have_collided.value) + # collision_interaction!(v_ode, u_ode, semi, semi.systems_have_collided.value) end - return dv_ode + return v_ode end # FluidSystems don't collide with each other -@inline function collision_interact!(dv_ode, v_ode, u_ode, system::FluidSystem, +@inline function collision_interact!(v_ode, u_ode, system::FluidSystem, neighbor::FluidSystem, semi; timer_str="") - return dv_ode + return v_ode end # BoundarySPHSystems don't collide with each other use RigidSPHSystem -@inline function collision_interact!(dv_ode, v_ode, u_ode, +@inline function collision_interact!(v_ode, u_ode, system::BoundarySPHSystem, neighbor::BoundarySPHSystem, semi; timer_str="") - return dv_ode + return v_ode end # @inline function set_collision_status!(semi, system) @@ -671,9 +667,8 @@ end # One can benchmark, e.g. the fluid-fluid interaction, with: # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; # @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); -@inline function collision_interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; - timer_str="") - dv = wrap_v(dv_ode, system, semi) +@inline function collision_interact!(v_ode, u_ode, system, neighbor, semi; timer_str="") + #dv = wrap_v(dv_ode, system, semi) # dv = wrap_u(dv_ode, system, semi) v_system = wrap_v(v_ode, system, semi) @@ -682,24 +677,25 @@ end v_neighbor = wrap_v(v_ode, neighbor, semi) u_neighbor = wrap_u(u_ode, neighbor, semi) nhs = get_neighborhood_search(system, neighbor, semi) + nhs_bnd = get_neighborhood_search(neighbor, neighbor, semi) @trixi_timeit timer() timer_str begin - collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, nhs, system, + collision_interact!(v_system, u_system, v_neighbor, u_neighbor, nhs, nhs_bnd, system, neighbor) end # set_collision_status!(semi, system) end -@inline function collision_interact!(dv, v_system, u_system, v_neighbor, u_neighbor, - nhs, system, neighbor) - return dv +@inline function collision_interact!(v_system, u_system, v_neighbor, u_neighbor, + nhs, nhs_bnd, system, neighbor) + return v_system end # Systems representing solids colliding with a boundary under the assumption the boundary is not influenced by the collision -@inline function collision_interact!(dv, v_particle_system, u_particle_system, +@inline function collision_interact!(v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, - neighborhood_search, + neighborhood_search, nhs_bnd, particle_system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}, neighbor_system::BoundarySPHSystem) @@ -709,8 +705,6 @@ end (; particle_spacing) = particle_system particle_radius = 0.5 * particle_spacing - collision_distance = neighbor_radius + particle_radius - system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -718,23 +712,87 @@ end particle_system.collision_u .= zeros(ndims(particle_system)) max_overlap = 0.0 collision_normal = zeros(ndims(particle_system)) + contact_point = zeros(ndims(particle_system)) + largest_overlap_coord = zeros(ndims(particle_system)) + + body_velocity = current_velocity(v_particle_system, particle_system, 1) # Detect the maximum overlap and calculate the normal vector at maximum overlap for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - overlap = collision_distance - distance - if overlap > max_overlap - max_overlap = overlap - # Normal vector pointing from particle to neighbor - collision_normal = pos_diff / distance + neighborhood_search, + parallel=false) do particle, neighbor, pos_diff, distance + + # Check axis-aligned bounding box overlap i.e. we assume the particles are actually the cell center of a cubic mesh + if all(abs(pos_diff[i]) <= (neighbor_radius + particle_radius) + for i in 1:ndims(particle_system)) + + bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) + overlap_coord = extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) + + # Calculate overlap on each axis and find the minimum to determine the collision normal + overlaps=[] + if abs(bnd_coord[2] - overlap_coord[2]) < neighbor_radius + particle_radius && bnd_coord[2] > 0.0 + overlaps = [abs(overlap_coord[1] - bnd_coord[1]), 0.0] + else + overlaps = [0.0, abs(overlap_coord[2] - bnd_coord[2])] + end + + overlaps = [(particle_radius + neighbor_radius) - abs(pos_diff[i]) + for i in 1:ndims(particle_system)] + # min_overlap, axis = findmin(overlaps) + # if min_overlap > max_overlap + # # max_overlap = min_overlap + # # collision_normal[axis] = sign(pos_diff[axis]) + # # bnd_coord = extract_svector(neighbor_system_coords, + # # Val(ndims(neighborhood_search)), neighbor) + # # largest_overlap_coord .= extract_svector(system_coords, + # # Val(ndims(neighborhood_search)), + # # particle) + # # contact_point = largest_overlap_coord + collision_normal * max_overlap + # normal_accumulator = zeros(Float64, ndims(particle_system)) + # bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) + + # # Calculate the surface normal using neighbors of the most overlapping boundary particle + # for bnd_neighbor in eachneighbor(bnd_coord, nhs_bnd) + # n_pos_diff = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), bnd_neighbor) - bnd_coord + # n_distance = norm(n_pos_diff) + # if n_distance > sqrt(eps()) + # normal_accumulator += normalize(n_pos_diff) * smoothing_kernel(neighbor_system, n_distance) + # end + # # println("n_pos_diff", n_pos_diff, " ", normalize(n_pos_diff)) + + # # println("normal_accumulator", normal_accumulator) + # end + + # collision_normal = normalize(normal_accumulator) + + if abs(overlaps[2]) > sqrt(eps()) + if overlaps[2] > max_overlap + collision_normal = [0.0, -1.0] + max_overlap = overlaps[2] + end + elseif abs(overlaps[1]) > sqrt(eps()) + if overlaps[1] > max_overlap + collision_normal = [-1.0, 0.0] + max_overlap = overlaps[1] + end + end + + # max_overlap = min_overlap + # overall_normal = normalize(pos_diff) # Calculate a normalized overall normal + # collision_normal = sign.(pos_diff) .* (overlaps ./ min_overlap) # Weighted direction for minimum overlap + # bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) + largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) + contact_point = largest_overlap_coord - 100.0 * collision_normal * max_overlap + end end if max_overlap > 0 particle_system.has_collided.value = true - initial_velocity_normal = dot(extract_svector(v_particle_system, particle_system, 1), collision_normal) + initial_velocity_normal = dot(body_velocity, collision_normal) # Coefficient of restitution (0 < e <= 1) # Tungsten Carbide = 0.7 to 1.0 @@ -743,30 +801,68 @@ end # Low Grade aluminium unhardened = 0.3 to 0.4 # Pure aluminium = 0.08 to 0.12 restitution_coefficient = 0.1 - if norm(initial_velocity_normal) < 10*sqrt(eps()) + if norm(initial_velocity_normal) < 10 * sqrt(eps()) restitution_coefficient = 0.0 end # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) - velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * collision_normal - - # we need to move in the opposite direction the same amount so times 2 - positional_correction = (2.0 * max_overlap + sqrt(eps())) * collision_normal - # if we don't move back far enough the current version gets unstable - movement_scale = max(1.0, 0.25 * particle_radius / norm(positional_correction)) - positional_correction_scaled = positional_correction * movement_scale - - particle_system.collision_u .= positional_correction_scaled - particle_system.collision_impulse .= velocity_change_normal + velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * + collision_normal + + println(particle_system) + println("initial_velocity_normal ", initial_velocity_normal) + println("collision_normal ", collision_normal) + + post_collision_velocity = body_velocity + velocity_change_normal + + # velocity_change = norm(v_particle_system[:, 1]) - norm(v_particle_system[:, 1] + velocity_change_normal) + # if abs(velocity_change) > 1 + # println("before: ", v_particle_system[:, 1]) + # println("after: ", v_particle_system[:, 1] + velocity_change_normal) + # println("imp: ", velocity_change_normal) + + # exit(-1) + # end + + # Calculate back movement time based on post-collision velocity + back_movement_time = (max_overlap + sqrt(eps())) / norm(body_velocity) + back_movement = post_collision_velocity * back_movement_time + + # positional_correction = max_overlap * collision_normal + + # # we need to move in the opposite direction the same amount so times 2 + # positional_correction = (2.0 * max_overlap + sqrt(eps())) * collision_normal + # # if we don't move back far enough the current version gets unstable + # movement_scale = max(1.0, 0.25 * particle_radius / norm(positional_correction)) + # positional_correction_scaled = positional_correction * movement_scale + + # position_correction = (contact_point - largest_overlap_coord) + # println("position_correction", position_correction) + # println("largest_overlap_coord", largest_overlap_coord) + println("contact_point", contact_point) + println("post_collision_velocity", post_collision_velocity) + # println("body_velocity", body_velocity) + # exit(-1) + + # if norm(position_correction) > 0.1 + # exit(-1) + # end + + position_correction = (contact_point - largest_overlap_coord) + + back_movement + # collision_impulse = velocity_change_normal + println("position_correction", position_correction) + println("back_movement", back_movement) # Apply the change uniformly across all particles for particle in each_moving_particle(particle_system) - current_normal_component = dot(dv[:, particle], collision_normal) * collision_normal - dv[:, particle] += -current_normal_component + v_particle_system[:, particle] = post_collision_velocity + # u_particle_system[:, particle] += position_correction + u_particle_system[:, particle] += [0.0, 0.5] end end - return dv + return v_particle_system end # collision influences both systems diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl index 9ff0451fd..85826aea4 100644 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -107,10 +107,6 @@ end # end @inline function current_velocity(v, system::RigidSPHSystem, particle) - if particle > n_moving_particles(system) - return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) - end - return extract_svector(v, system, particle) end From 47c570f9542b0ddb8cdfaf323e8bad9163cdd329 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 17 Apr 2024 18:36:41 +0200 Subject: [PATCH 08/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 48 +++++++++++-- src/general/semidiscretization.jl | 71 +++++++++++-------- .../dummy_particles/dummy_particles.jl | 10 ++- src/schemes/solid/rigid_body_sph/system.jl | 6 ++ 4 files changed, 96 insertions(+), 39 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index 2c68ce3e4..b45ca0ca6 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -14,7 +14,7 @@ spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -tspan = (0.0, 0.4) +tspan = (0.0, 2.0) # Boundary geometry and initial fluid particle positions initial_fluid_size = (2.0, 0.5) @@ -104,19 +104,18 @@ solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(boundary_system, solid_system_2) +semi = Semidiscretization(boundary_system, solid_system_2, fluid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", write_meta_data=true) -callbacks = CallbackSet(info_callback, saving_callback) function collision(vu, integrator, semi, t) v_ode, u_ode = vu.x - println("vollision at ", t) + # println("vollision at ", t) TrixiParticles.collision_interaction!(v_ode, u_ode, semi) @@ -141,8 +140,43 @@ function collision(vu, integrator, semi, t) # end end +function collision2(vu, integrator, semi, t) + v_ode, u_ode = vu.x + + # println("vollision at ", t) + + TrixiParticles.collision_interaction!(v_ode, u_ode, semi, u_only=true) + + # TrixiParticles.foreach_system(semi) do system + # v = TrixiParticles.wrap_v(v_ode, system, semi) + # u = TrixiParticles.wrap_u(u_ode, system, semi) + + # if system isa RigidSPHSystem && system.has_collided.value + # velocity_change = norm(v[:, 1]) - norm(v[:, 1] + system.collision_impulse) + # if abs(velocity_change) > 1 + # println("before: ", v[:, 1]) + # println("after: ", v[:, 1] + system.collision_impulse) + # println("imp: ", system.collision_impulse) + + # exit(-1) + # end + # for particle in TrixiParticles.each_moving_particle(system) + # v[:, particle] += system.collision_impulse + # u[:, particle] += system.collision_u + # end + # end + # end +end + # Use a Runge-Kutta method with automatic (error based) time step size control. -sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), - abstol=1e-7, # Default abstol is 1e-6 - reltol=1e-5, # Default reltol is 1e-3 +# sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), +# abstol=1e-7, # Default abstol is 1e-6 +# reltol=1e-5, # Default reltol is 1e-3 +# save_everystep=false, callback=callbacks); +stepsize_callback = StepsizeCallback(cfl=0.5) + +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false;step_limiter! =collision), + dt=1.0, # This is overwritten by the stepsize callback save_everystep=false, callback=callbacks); diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index b236c7b95..b6afbd28c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -368,14 +368,16 @@ function drift!(du_ode, v_ode, u_ode, semi, t) foreach_system(semi) do system du = wrap_u(du_ode, system, semi) v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) @threaded for particle in each_moving_particle(system) # This can be dispatched per system add_velocity!(du, v, particle, system) end - if system isa RigidSPHSystem && system.has_collided.value - println("body speed", v[:,1]) - end + # if system isa RigidSPHSystem + # println("body speed", v[:,1], " pos ", u[:,1]) + # # println("in drift:", u[:,1]) + # end end end @@ -397,7 +399,6 @@ end @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du function kick!(dv_ode, v_ode, u_ode, semi, t) - println("kick at ", t) @trixi_timeit timer() "kick!" begin @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode) @@ -593,24 +594,24 @@ end end # start of recursion -function collision_interaction!(v_ode, u_ode, semi) - collision_interaction!(v_ode, u_ode, semi, true) +function collision_interaction!(v_ode, u_ode, semi; u_only=false) + collision_interaction!(v_ode, u_ode, semi, true, u_only) return v_ode end -# @inline function reset_system_collision_status(semi) -# foreach_system(semi) do system -# if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem -# system.has_collided.value = false -# end -# end -# end +@inline function reset_system_collision_status(semi) + foreach_system(semi) do system + if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem + system.has_collided.value = false + end + end +end -function collision_interaction!(v_ode, u_ode, semi, systems_have_collided) +function collision_interaction!(v_ode, u_ode, semi, systems_have_collided, u_only) if systems_have_collided semi.systems_have_collided.value = false - # reset_system_collision_status(semi) + reset_system_collision_status(semi) # Call `interact!` for each pair of systems foreach_system(semi) do system @@ -628,7 +629,7 @@ function collision_interaction!(v_ode, u_ode, semi, systems_have_collided) timer_str = "" end - collision_interact!(v_ode, u_ode, system, neighbor, semi, + collision_interact!(v_ode, u_ode, system, neighbor, semi, u_only, timer_str=timer_str) end end @@ -643,14 +644,14 @@ end # FluidSystems don't collide with each other @inline function collision_interact!(v_ode, u_ode, system::FluidSystem, - neighbor::FluidSystem, semi; timer_str="") + neighbor::FluidSystem, semi, u_only; timer_str="") return v_ode end # BoundarySPHSystems don't collide with each other use RigidSPHSystem @inline function collision_interact!(v_ode, u_ode, system::BoundarySPHSystem, neighbor::BoundarySPHSystem, - semi; timer_str="") + semi, u_only; timer_str="") return v_ode end @@ -667,7 +668,7 @@ end # One can benchmark, e.g. the fluid-fluid interaction, with: # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; # @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); -@inline function collision_interact!(v_ode, u_ode, system, neighbor, semi; timer_str="") +@inline function collision_interact!(v_ode, u_ode, system, neighbor, semi, u_only; timer_str="") #dv = wrap_v(dv_ode, system, semi) # dv = wrap_u(dv_ode, system, semi) @@ -681,14 +682,14 @@ end @trixi_timeit timer() timer_str begin collision_interact!(v_system, u_system, v_neighbor, u_neighbor, nhs, nhs_bnd, system, - neighbor) + neighbor, u_only) end # set_collision_status!(semi, system) end @inline function collision_interact!(v_system, u_system, v_neighbor, u_neighbor, - nhs, nhs_bnd, system, neighbor) + nhs, nhs_bnd, system, neighbor, u_only) return v_system end @@ -698,7 +699,7 @@ end neighborhood_search, nhs_bnd, particle_system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}, - neighbor_system::BoundarySPHSystem) + neighbor_system::BoundarySPHSystem, u_only) (; particle_spacing) = neighbor_system neighbor_radius = 0.5 * particle_spacing @@ -727,8 +728,11 @@ end if all(abs(pos_diff[i]) <= (neighbor_radius + particle_radius) for i in 1:ndims(particle_system)) + # println("pos_diff", pos_diff) + bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) overlap_coord = extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) + # println("coord", bnd_coord, " ", overlap_coord) # Calculate overlap on each axis and find the minimum to determine the collision normal overlaps=[] @@ -771,11 +775,15 @@ end if overlaps[2] > max_overlap collision_normal = [0.0, -1.0] max_overlap = overlaps[2] + largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) + contact_point = largest_overlap_coord - collision_normal * max_overlap end elseif abs(overlaps[1]) > sqrt(eps()) if overlaps[1] > max_overlap collision_normal = [-1.0, 0.0] max_overlap = overlaps[1] + largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) + contact_point = largest_overlap_coord - collision_normal * max_overlap end end @@ -783,8 +791,7 @@ end # overall_normal = normalize(pos_diff) # Calculate a normalized overall normal # collision_normal = sign.(pos_diff) .* (overlaps ./ min_overlap) # Weighted direction for minimum overlap # bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) - contact_point = largest_overlap_coord - 100.0 * collision_normal * max_overlap + end end @@ -800,8 +807,10 @@ end # High Grade aluminium hardened = 0.7 to 0.8 # Low Grade aluminium unhardened = 0.3 to 0.4 # Pure aluminium = 0.08 to 0.12 - restitution_coefficient = 0.1 - if norm(initial_velocity_normal) < 10 * sqrt(eps()) + restitution_coefficient = 0.5 + println("initial_velocity_normal", initial_velocity_normal, " max ", max_overlap) + # todo: this actually depends on the timestep + if norm(initial_velocity_normal) < 1e-2 || max_overlap < 1e-6 restitution_coefficient = 0.0 end @@ -809,8 +818,6 @@ end velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * collision_normal - println(particle_system) - println("initial_velocity_normal ", initial_velocity_normal) println("collision_normal ", collision_normal) post_collision_velocity = body_velocity + velocity_change_normal @@ -857,9 +864,13 @@ end # Apply the change uniformly across all particles for particle in each_moving_particle(particle_system) v_particle_system[:, particle] = post_collision_velocity - # u_particle_system[:, particle] += position_correction - u_particle_system[:, particle] += [0.0, 0.5] + u_particle_system[:, particle] .+= position_correction + # if min_y > u_particle_system[2, particle] + # min_y = u_particle_system[2, particle] + # end + # u_particle_system[:, particle] += [0.0, 0.25*neighbor_radius] end + # println("miny ", min_y) end return v_particle_system diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 27c7f924d..e0961f7e8 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -378,12 +378,18 @@ end particles=eachparticle(system)) do particle, neighbor, pos_diff, distance density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) + kernel_weight = smoothing_kernel(boundary_model, distance) + # if system isa RigidSPHSystem + # resulting_acc = neighbor_system.acceleration - 100.0 * system.acceleration + # else resulting_acc = neighbor_system.acceleration - current_acceleration(system, particle) + # end - kernel_weight = smoothing_kernel(boundary_model, distance) - + if system isa RigidSPHSystem + pressure[particle] +=100000*kernel_weight + end pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, neighbor) + dot(resulting_acc, density_neighbor * pos_diff)) * diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl index 85826aea4..44de55356 100644 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ b/src/schemes/solid/rigid_body_sph/system.jl @@ -197,3 +197,9 @@ end function viscosity_model(system::RigidSPHSystem) return system.boundary_model.viscosity end + +function calculate_dt(v_ode, u_ode, cfl_number, system::RigidSPHSystem) + (; acceleration, particle_spacing) = system + + return cfl_number * particle_spacing/maximum(abs.(acceleration)) +end From ba81acbf62361ead2d4285671f3452e803e98a57 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 18 Apr 2024 01:43:39 +0200 Subject: [PATCH 09/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 36 +++---------------- src/general/semidiscretization.jl | 12 +++---- .../dummy_particles/dummy_particles.jl | 14 +++++--- 3 files changed, 20 insertions(+), 42 deletions(-) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index b45ca0ca6..7cff6f3e0 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -4,7 +4,7 @@ using LinearAlgebra # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.015 +fluid_particle_spacing = 0.01 solid_particle_spacing = fluid_particle_spacing # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model @@ -14,7 +14,7 @@ spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -tspan = (0.0, 2.0) +tspan = (0.0, 4.0) # Boundary geometry and initial fluid particle positions initial_fluid_size = (2.0, 0.5) @@ -39,7 +39,7 @@ sphere1_density = 600.0 sphere2_density = 3000 sphere1_center = (0.5, 1.1) -sphere2_center = (1.5, 0.8) +sphere2_center = (1.0, 0.8) sphere1 = SphereShape(solid_particle_spacing, sphere1_radius, sphere1_center, sphere1_density, sphere_type=VoxelSphere()) sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, @@ -47,7 +47,7 @@ sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, # ========================================================================================== # ==== Fluid -fluid_smoothing_length = 3.0 * fluid_particle_spacing +fluid_smoothing_length = 3.5 * fluid_particle_spacing fluid_smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() @@ -140,34 +140,6 @@ function collision(vu, integrator, semi, t) # end end -function collision2(vu, integrator, semi, t) - v_ode, u_ode = vu.x - - # println("vollision at ", t) - - TrixiParticles.collision_interaction!(v_ode, u_ode, semi, u_only=true) - - # TrixiParticles.foreach_system(semi) do system - # v = TrixiParticles.wrap_v(v_ode, system, semi) - # u = TrixiParticles.wrap_u(u_ode, system, semi) - - # if system isa RigidSPHSystem && system.has_collided.value - # velocity_change = norm(v[:, 1]) - norm(v[:, 1] + system.collision_impulse) - # if abs(velocity_change) > 1 - # println("before: ", v[:, 1]) - # println("after: ", v[:, 1] + system.collision_impulse) - # println("imp: ", system.collision_impulse) - - # exit(-1) - # end - # for particle in TrixiParticles.each_moving_particle(system) - # v[:, particle] += system.collision_impulse - # u[:, particle] += system.collision_u - # end - # end - # end -end - # Use a Runge-Kutta method with automatic (error based) time step size control. # sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), # abstol=1e-7, # Default abstol is 1e-6 diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index b6afbd28c..8bbfff56f 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -808,7 +808,7 @@ end # Low Grade aluminium unhardened = 0.3 to 0.4 # Pure aluminium = 0.08 to 0.12 restitution_coefficient = 0.5 - println("initial_velocity_normal", initial_velocity_normal, " max ", max_overlap) + # println("initial_velocity_normal", initial_velocity_normal, " max ", max_overlap) # todo: this actually depends on the timestep if norm(initial_velocity_normal) < 1e-2 || max_overlap < 1e-6 restitution_coefficient = 0.0 @@ -818,7 +818,7 @@ end velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * collision_normal - println("collision_normal ", collision_normal) + # println("collision_normal ", collision_normal) post_collision_velocity = body_velocity + velocity_change_normal @@ -846,8 +846,8 @@ end # position_correction = (contact_point - largest_overlap_coord) # println("position_correction", position_correction) # println("largest_overlap_coord", largest_overlap_coord) - println("contact_point", contact_point) - println("post_collision_velocity", post_collision_velocity) + # println("contact_point", contact_point) + # println("post_collision_velocity", post_collision_velocity) # println("body_velocity", body_velocity) # exit(-1) @@ -858,8 +858,8 @@ end position_correction = (contact_point - largest_overlap_coord) + back_movement # collision_impulse = velocity_change_normal - println("position_correction", position_correction) - println("back_movement", back_movement) + # println("position_correction", position_correction) + # println("back_movement", back_movement) # Apply the change uniformly across all particles for particle in each_moving_particle(particle_system) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index e0961f7e8..dfb642eb0 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -329,7 +329,7 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, + system_coords, neighbor_coords, v, v_neighbor_system, nhs) end @@ -367,7 +367,7 @@ end @inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system::FluidSystem, - system_coords, neighbor_coords, + system_coords, neighbor_coords, v_particle_system, v_neighbor_system, neighborhood_search) (; pressure, cache, viscosity) = boundary_model @@ -388,7 +388,12 @@ end # end if system isa RigidSPHSystem - pressure[particle] +=100000*kernel_weight + # pressure[particle] += 4000 * kernel_weight + pressure[particle] += (10000 + 0.5 * density_neighbor * norm(v_particle_system[1:ndims(system),1])^2) * kernel_weight + # println("old ", 0.5* density_neighbor *norm(v_particle_system[:,1])^2 * kernel_weight, " ", 4000 * kernel_weight) + # println("new ", 0.5* density_neighbor *norm(v_particle_system[1:ndims(system),1]- v_neighbor_system[1:ndims(system),neighbor])^2 * kernel_weight, " ", 4000 * kernel_weight) + # else + # pressure[particle] += (0.5 * density_neighbor * norm(v_neighbor_system[1:ndims(system),neighbor])^2) * kernel_weight end pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, neighbor) + @@ -397,6 +402,7 @@ end cache.volume[particle] += kernel_weight + # why is this done here? compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) end @@ -409,7 +415,7 @@ end end @inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, + system_coords, neighbor_coords, v_particle_system, v_neighbor_system, neighborhood_search) return boundary_model end From 72f9ad27e94ad97a2f626241e1334dbe3c4317a1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 22 Apr 2024 14:28:22 +0200 Subject: [PATCH 10/10] update --- examples/fsi/falling_rigid_spheres_2d.jl | 154 ------- src/TrixiParticles.jl | 2 +- src/callbacks/post_process.jl | 8 +- src/callbacks/solution_saving.jl | 67 ++- src/general/general.jl | 4 - src/general/semidiscretization.jl | 351 ++------------ .../dummy_particles/dummy_particles.jl | 123 +++-- src/schemes/boundary/rhs.jl | 6 +- src/schemes/boundary/system.jl | 8 +- .../fluid/weakly_compressible_sph/rhs.jl | 37 -- src/schemes/schemes.jl | 6 - src/schemes/solid/rigid_body_sph/rhs.jl | 118 ----- .../solid/rigid_body_sph/rigid_body_sph.jl | 1 - src/schemes/solid/rigid_body_sph/system.jl | 205 --------- src/visualization/write2vtk.jl | 7 - test/callbacks/callbacks.jl | 1 + test/callbacks/solution_saving.jl | 68 +++ test/validation/validation.jl | 10 +- .../dam_break_2d/validation_dam_break_2d.jl | 8 +- .../validation_reference_edac_0015.json | 412 ++++++++--------- .../validation_reference_wcsph_0015.json | 428 +++++++++--------- validation/validation_util.jl | 46 ++ 22 files changed, 722 insertions(+), 1348 deletions(-) delete mode 100644 examples/fsi/falling_rigid_spheres_2d.jl delete mode 100644 src/schemes/solid/rigid_body_sph/rhs.jl delete mode 100644 src/schemes/solid/rigid_body_sph/rigid_body_sph.jl delete mode 100644 src/schemes/solid/rigid_body_sph/system.jl create mode 100644 test/callbacks/solution_saving.jl diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl deleted file mode 100644 index 7cff6f3e0..000000000 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ /dev/null @@ -1,154 +0,0 @@ -using TrixiParticles -using OrdinaryDiffEq -using LinearAlgebra - -# ========================================================================================== -# ==== Resolution -fluid_particle_spacing = 0.01 -solid_particle_spacing = fluid_particle_spacing - -# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model -boundary_layers = 3 -spacing_ratio = 1 - -# ========================================================================================== -# ==== Experiment Setup -gravity = 9.81 -tspan = (0.0, 4.0) - -# Boundary geometry and initial fluid particle positions -initial_fluid_size = (2.0, 0.5) -tank_size = (2.0, 2.0) - -fluid_density = 1000.0 -sound_speed = 150 -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=1) - -tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - n_layers=boundary_layers, spacing_ratio=spacing_ratio, - faces=(true, true, true, false), - acceleration=(0.0, -gravity), state_equation=state_equation) - -sphere1_radius = 0.3 -sphere2_radius = 0.2 -# wood -sphere1_density = 600.0 -# steel -#sphere2_density = 7700 -sphere2_density = 3000 - -sphere1_center = (0.5, 1.1) -sphere2_center = (1.0, 0.8) -sphere1 = SphereShape(solid_particle_spacing, sphere1_radius, sphere1_center, - sphere1_density, sphere_type=VoxelSphere()) -sphere2 = SphereShape(solid_particle_spacing, sphere2_radius, sphere2_center, - sphere2_density, sphere_type=VoxelSphere()) - -# ========================================================================================== -# ==== Fluid -fluid_smoothing_length = 3.5 * fluid_particle_spacing -fluid_smoothing_kernel = WendlandC2Kernel{2}() - -fluid_density_calculator = ContinuityDensity() -viscosity = ViscosityAdami(nu=1e-4) -density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) - -fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, - state_equation, fluid_smoothing_kernel, - fluid_smoothing_length, viscosity=viscosity, - density_diffusion=density_diffusion, - acceleration=(0.0, -gravity)) - -# ========================================================================================== -# ==== Boundary -boundary_density_calculator = AdamiPressureExtrapolation() -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - state_equation=state_equation, - boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length, viscosity=viscosity) - -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, - particle_spacing=fluid_particle_spacing) - -# ========================================================================================== -# ==== Solid - -# For the FSI we need the hydrodynamic masses and densities in the solid boundary model -hydrodynamic_densites_1 = fluid_density * ones(size(sphere1.density)) -hydrodynamic_masses_1 = hydrodynamic_densites_1 * solid_particle_spacing^ndims(fluid_system) - -solid_boundary_model_1 = BoundaryModelDummyParticles(hydrodynamic_densites_1, - hydrodynamic_masses_1, - state_equation=state_equation, - boundary_density_calculator, - fluid_smoothing_kernel, - fluid_smoothing_length, viscosity=viscosity) - -hydrodynamic_densites_2 = fluid_density * ones(size(sphere2.density)) -hydrodynamic_masses_2 = hydrodynamic_densites_2 * solid_particle_spacing^ndims(fluid_system) - -solid_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2, - hydrodynamic_masses_2, - state_equation=state_equation, - boundary_density_calculator, - fluid_smoothing_kernel, - fluid_smoothing_length,viscosity=viscosity) - -solid_system_1 = RigidSPHSystem(sphere1; boundary_model=solid_boundary_model_1, - acceleration=(0.0, -gravity), - particle_spacing=fluid_particle_spacing) -solid_system_2 = RigidSPHSystem(sphere2; boundary_model=solid_boundary_model_2, - acceleration=(0.0, -gravity), - particle_spacing=fluid_particle_spacing) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(boundary_system, solid_system_2, fluid_system) -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=10) -saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", - write_meta_data=true) - - -function collision(vu, integrator, semi, t) - v_ode, u_ode = vu.x - - # println("vollision at ", t) - - TrixiParticles.collision_interaction!(v_ode, u_ode, semi) - - # TrixiParticles.foreach_system(semi) do system - # v = TrixiParticles.wrap_v(v_ode, system, semi) - # u = TrixiParticles.wrap_u(u_ode, system, semi) - - # if system isa RigidSPHSystem && system.has_collided.value - # velocity_change = norm(v[:, 1]) - norm(v[:, 1] + system.collision_impulse) - # if abs(velocity_change) > 1 - # println("before: ", v[:, 1]) - # println("after: ", v[:, 1] + system.collision_impulse) - # println("imp: ", system.collision_impulse) - - # exit(-1) - # end - # for particle in TrixiParticles.each_moving_particle(system) - # v[:, particle] += system.collision_impulse - # u[:, particle] += system.collision_u - # end - # end - # end -end - -# Use a Runge-Kutta method with automatic (error based) time step size control. -# sol = solve(ode, RDPK3SpFSAL49(;stage_limiter! =collision), -# abstol=1e-7, # Default abstol is 1e-6 -# reltol=1e-5, # Default reltol is 1e-3 -# save_everystep=false, callback=callbacks); -stepsize_callback = StepsizeCallback(cfl=0.5) - -callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback) - -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false;step_limiter! =collision), - dt=1.0, # This is overwritten by the stepsize callback - save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 2298cc423..b7c0c5fe7 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -42,7 +42,7 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, RigidSPHSystem + BoundarySPHSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback export ContinuityDensity, SummationDensity diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 243789433..25a8d051b 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -96,12 +96,12 @@ function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=tru exclude_boundary, funcs, filename, output_directory, append_timestamp, write_csv, write_json) if dt > 0 - # Add a `tstop` every `dt`, and save the final solution. + # Add a `tstop` every `dt`, and save the final solution return PeriodicCallback(post_callback, dt, initialize=initialize_postprocess_callback!, save_positions=(false, false), final_affect=true) else - # The first one is the condition, the second the affect! + # The first one is the `condition`, the second the `affect!` return DiscreteCallback(post_callback, post_callback, save_positions=(false, false), initialize=initialize_postprocess_callback!) @@ -214,7 +214,7 @@ function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integra return nothing end -# condition with interval +# `condition` with interval function (pp::PostprocessCallback)(u, t, integrator) (; interval) = pp @@ -260,7 +260,7 @@ function (pp::PostprocessCallback)(integrator) write_postprocess_callback(pp) end - # Tell OrdinaryDiffEq that u has not been modified + # Tell OrdinaryDiffEq that `u` has not been modified u_modified!(integrator, false) end diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 119f80e80..b83e9d4ef 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -104,15 +104,16 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0, -1) if length(save_times) > 0 - return PresetTimeCallback(save_times, solution_callback) + # See the large comment below for an explanation why we use `finalize` here + return PresetTimeCallback(save_times, solution_callback, finalize=solution_callback) elseif dt > 0 - # Add a `tstop` every `dt`, and save the final solution. + # Add a `tstop` every `dt`, and save the final solution return PeriodicCallback(solution_callback, dt, initialize=initialize_save_cb!, save_positions=(false, false), final_affect=save_final_solution) else - # The first one is the condition, the second the affect! + # The first one is the `condition`, the second the `affect!` return DiscreteCallback(solution_callback, solution_callback, save_positions=(false, false), initialize=initialize_save_cb!) @@ -132,19 +133,19 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in # Save initial solution if solution_callback.save_initial_solution - # Update systems to compute quantities like density and pressure. + # Update systems to compute quantities like density and pressure semi = integrator.p v_ode, u_ode = u.x update_systems_and_nhs(v_ode, u_ode, semi, t) - # Apply the callback. + # Apply the callback solution_callback(integrator) end return nothing end -# condition +# `condition` function (solution_callback::SolutionSavingCallback)(u, t, integrator) (; interval, save_final_solution) = solution_callback @@ -152,7 +153,7 @@ function (solution_callback::SolutionSavingCallback)(u, t, integrator) save_final_solution=save_final_solution) end -# affect! +# `affect!` function (solution_callback::SolutionSavingCallback)(integrator) (; interval, output_directory, custom_quantities, write_meta_data, verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback @@ -183,12 +184,30 @@ function (solution_callback::SolutionSavingCallback)(integrator) max_coordinates=max_coordinates, custom_quantities...) - # Tell OrdinaryDiffEq that u has not been modified + # Tell OrdinaryDiffEq that `u` has not been modified u_modified!(integrator, false) return nothing end +# `finalize` +# This is a hack to make dispatch on a `PresetTimeCallback` possible. +# +# The type of the returned `DiscreteCallback` is +# `DiscreteCallback{typeof(condition), typeof(affect!), typeof(initialize), typeof(finalize)}`. +# For the `PeriodicCallback`, `typeof(affect!)` contains the type of the +# `SolutionSavingCallback`. The `PresetTimeCallback` uses anonymous functions as `condition` +# and `affect!`, so this doesn't work here. +# +# This hacky workaround makes use of the type parameter `typeof(finalize)` above. +# It's set to `FINALIZE_DEFAULT` by default in the `PresetTimeCallback`, which is a function +# that just returns `nothing`. +# Instead, we pass the `SolutionSavingCallback` as `finalize`, and define it to also just +# return `nothing` when called as `initialize`. +function (finalize::SolutionSavingCallback)(c, u, t, integrator) + return nothing +end + function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time @@ -205,6 +224,14 @@ function Base.show(io::IO, print(io, "SolutionSavingCallback(dt=", solution_saving.interval, ")") end +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, <:Any, <:Any, <:SolutionSavingCallback}) + @nospecialize cb # reduce precompilation time + + solution_saving = cb.finalize + print(io, "SolutionSavingCallback(save_times=", solution_saving.save_times, ")") +end + function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time @@ -229,6 +256,30 @@ function Base.show(io::IO, ::MIME"text/plain", end end +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:Any, <:Any, <:SolutionSavingCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + solution_saving = cb.finalize + cq = collect(solution_saving.custom_quantities) + + setup = [ + "save_times" => solution_saving.save_times, + "custom quantities" => isempty(cq) ? nothing : cq, + "save initial solution" => solution_saving.save_initial_solution ? + "yes" : "no", + "save final solution" => solution_saving.save_final_solution ? "yes" : + "no", + "output directory" => abspath(solution_saving.output_directory), + "prefix" => solution_saving.prefix, + ] + summary_box(io, "SolutionSavingCallback", setup) + end +end + function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:PeriodicCallbackAffect{<:SolutionSavingCallback}}) diff --git a/src/general/general.jl b/src/general/general.jl index 9576cf2c7..6651c03c5 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -15,10 +15,6 @@ timer_name(::BoundarySystem) = "boundary" return du end -mutable struct MutableBool - value::Bool -end - # Note that `semidiscretization.jl` depends on the system types and has to be # included later. # `density_calculators.jl` needs to be included before `corrections.jl`. diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 8bbfff56f..ada9275e4 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -49,7 +49,6 @@ struct Semidiscretization{S, RU, RV, NS} ranges_u :: RU ranges_v :: RV neighborhood_searches :: NS - systems_have_collided :: MutableBool # Dispatch at `systems` to distinguish this constructor from the one below when # 4 systems are passed. @@ -57,8 +56,7 @@ struct Semidiscretization{S, RU, RV, NS} function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches) new{typeof(systems), typeof(ranges_u), typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v, - neighborhood_searches, - MutableBool(false)) + neighborhood_searches) end end @@ -158,8 +156,7 @@ end return compact_support(smoothing_kernel, smoothing_length) end -@inline function compact_support(system::Union{TotalLagrangianSPHSystem, BoundarySPHSystem, - RigidSPHSystem}, +@inline function compact_support(system::Union{TotalLagrangianSPHSystem, BoundarySPHSystem}, neighbor) return compact_support(system, system.boundary_model, neighbor) end @@ -355,10 +352,6 @@ function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) end -function calculate_dt(v_ode, u_ode, cfl_number, system) - return Inf -end - function drift!(du_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "drift!" begin @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode) @@ -368,21 +361,13 @@ function drift!(du_ode, v_ode, u_ode, semi, t) foreach_system(semi) do system du = wrap_u(du_ode, system, semi) v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) @threaded for particle in each_moving_particle(system) # This can be dispatched per system add_velocity!(du, v, particle, system) end - # if system isa RigidSPHSystem - # println("body speed", v[:,1], " pos ", u[:,1]) - # # println("in drift:", u[:,1]) - # end end end - - # @trixi_timeit timer() "collision interaction" collision_interaction!(du_ode, v_ode, - # u_ode, semi) end return du_ode @@ -485,7 +470,6 @@ function add_source_terms!(dv_ode, v_ode, u_ode, semi) end @inline source_terms(system) = nothing -@inline source_terms(system::RigidSPHSystem) = nothing @inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms @inline add_acceleration!(dv, particle, system) = dv @@ -593,320 +577,39 @@ end end end -# start of recursion -function collision_interaction!(v_ode, u_ode, semi; u_only=false) - collision_interaction!(v_ode, u_ode, semi, true, u_only) - - return v_ode -end - -@inline function reset_system_collision_status(semi) - foreach_system(semi) do system - if system isa RigidSPHSystem || system isa TotalLagrangianSPHSystem - system.has_collided.value = false - end - end -end - -function collision_interaction!(v_ode, u_ode, semi, systems_have_collided, u_only) - if systems_have_collided - semi.systems_have_collided.value = false - reset_system_collision_status(semi) - - # Call `interact!` for each pair of systems - foreach_system(semi) do system - foreach_system(semi) do neighbor - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor, semi) - - # the same system does not collide with its self - if system_index != neighbor_index - # Construct string for the interactions timer. - # Avoid allocations from string construction when no timers are used. - if timeit_debug_enabled() - timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" - else - timer_str = "" - end - - collision_interact!(v_ode, u_ode, system, neighbor, semi, u_only, - timer_str=timer_str) - end - end - end - - # if systems have collided we need to call this again since there can be more collisions - # collision_interaction!(v_ode, u_ode, semi, semi.systems_have_collided.value) - end - - return v_ode -end - -# FluidSystems don't collide with each other -@inline function collision_interact!(v_ode, u_ode, system::FluidSystem, - neighbor::FluidSystem, semi, u_only; timer_str="") - return v_ode -end - -# BoundarySPHSystems don't collide with each other use RigidSPHSystem -@inline function collision_interact!(v_ode, u_ode, - system::BoundarySPHSystem, neighbor::BoundarySPHSystem, - semi, u_only; timer_str="") - return v_ode -end - -# @inline function set_collision_status!(semi, system) -# return semi -# end - -# @inline function set_collision_status!(semi, system::Union{RigidSPHSystem, TotalLagrangianSPHSystem}) -# semi.systems_have_collided.value = semi.systems_have_collided.value || system.has_collided.value -# return semi -# end - -# Function barrier to make benchmarking interactions easier. -# One can benchmark, e.g. the fluid-fluid interaction, with: -# dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; -# @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); -@inline function collision_interact!(v_ode, u_ode, system, neighbor, semi, u_only; timer_str="") - #dv = wrap_v(dv_ode, system, semi) - # dv = wrap_u(dv_ode, system, semi) - - v_system = wrap_v(v_ode, system, semi) - u_system = wrap_u(u_ode, system, semi) - - v_neighbor = wrap_v(v_ode, neighbor, semi) - u_neighbor = wrap_u(u_ode, neighbor, semi) - nhs = get_neighborhood_search(system, neighbor, semi) - nhs_bnd = get_neighborhood_search(neighbor, neighbor, semi) - - @trixi_timeit timer() timer_str begin - collision_interact!(v_system, u_system, v_neighbor, u_neighbor, nhs, nhs_bnd, system, - neighbor, u_only) - end - - # set_collision_status!(semi, system) +# NHS updates +function nhs_coords(system::FluidSystem, + neighbor::FluidSystem, u) + return current_coordinates(u, neighbor) end -@inline function collision_interact!(v_system, u_system, v_neighbor, u_neighbor, - nhs, nhs_bnd, system, neighbor, u_only) - return v_system +function nhs_coords(system::FluidSystem, + neighbor::TotalLagrangianSPHSystem, u) + return current_coordinates(u, neighbor) end -# Systems representing solids colliding with a boundary under the assumption the boundary is not influenced by the collision -@inline function collision_interact!(v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search, nhs_bnd, - particle_system::Union{RigidSPHSystem, - TotalLagrangianSPHSystem}, - neighbor_system::BoundarySPHSystem, u_only) - (; particle_spacing) = neighbor_system - neighbor_radius = 0.5 * particle_spacing - - (; particle_spacing) = particle_system - particle_radius = 0.5 * particle_spacing - - system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - - particle_system.collision_impulse .= zeros(ndims(particle_system)) - particle_system.collision_u .= zeros(ndims(particle_system)) - max_overlap = 0.0 - collision_normal = zeros(ndims(particle_system)) - contact_point = zeros(ndims(particle_system)) - largest_overlap_coord = zeros(ndims(particle_system)) - - body_velocity = current_velocity(v_particle_system, particle_system, 1) - - # Detect the maximum overlap and calculate the normal vector at maximum overlap - for_particle_neighbor(particle_system, neighbor_system, system_coords, - neighbor_system_coords, - neighborhood_search, - parallel=false) do particle, neighbor, pos_diff, distance - - # Check axis-aligned bounding box overlap i.e. we assume the particles are actually the cell center of a cubic mesh - if all(abs(pos_diff[i]) <= (neighbor_radius + particle_radius) - for i in 1:ndims(particle_system)) - - # println("pos_diff", pos_diff) - - bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - overlap_coord = extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) - # println("coord", bnd_coord, " ", overlap_coord) - - # Calculate overlap on each axis and find the minimum to determine the collision normal - overlaps=[] - if abs(bnd_coord[2] - overlap_coord[2]) < neighbor_radius + particle_radius && bnd_coord[2] > 0.0 - overlaps = [abs(overlap_coord[1] - bnd_coord[1]), 0.0] - else - overlaps = [0.0, abs(overlap_coord[2] - bnd_coord[2])] - end - - overlaps = [(particle_radius + neighbor_radius) - abs(pos_diff[i]) - for i in 1:ndims(particle_system)] - # min_overlap, axis = findmin(overlaps) - # if min_overlap > max_overlap - # # max_overlap = min_overlap - # # collision_normal[axis] = sign(pos_diff[axis]) - # # bnd_coord = extract_svector(neighbor_system_coords, - # # Val(ndims(neighborhood_search)), neighbor) - # # largest_overlap_coord .= extract_svector(system_coords, - # # Val(ndims(neighborhood_search)), - # # particle) - # # contact_point = largest_overlap_coord + collision_normal * max_overlap - # normal_accumulator = zeros(Float64, ndims(particle_system)) - # bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - - # # Calculate the surface normal using neighbors of the most overlapping boundary particle - # for bnd_neighbor in eachneighbor(bnd_coord, nhs_bnd) - # n_pos_diff = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), bnd_neighbor) - bnd_coord - # n_distance = norm(n_pos_diff) - # if n_distance > sqrt(eps()) - # normal_accumulator += normalize(n_pos_diff) * smoothing_kernel(neighbor_system, n_distance) - # end - # # println("n_pos_diff", n_pos_diff, " ", normalize(n_pos_diff)) - - # # println("normal_accumulator", normal_accumulator) - # end - - # collision_normal = normalize(normal_accumulator) - - if abs(overlaps[2]) > sqrt(eps()) - if overlaps[2] > max_overlap - collision_normal = [0.0, -1.0] - max_overlap = overlaps[2] - largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) - contact_point = largest_overlap_coord - collision_normal * max_overlap - end - elseif abs(overlaps[1]) > sqrt(eps()) - if overlaps[1] > max_overlap - collision_normal = [-1.0, 0.0] - max_overlap = overlaps[1] - largest_overlap_coord .= extract_svector(system_coords, Val(ndims(neighborhood_search)), particle) - contact_point = largest_overlap_coord - collision_normal * max_overlap - end - end - - # max_overlap = min_overlap - # overall_normal = normalize(pos_diff) # Calculate a normalized overall normal - # collision_normal = sign.(pos_diff) .* (overlaps ./ min_overlap) # Weighted direction for minimum overlap - # bnd_coord = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - - - end - end - - if max_overlap > 0 - particle_system.has_collided.value = true - - initial_velocity_normal = dot(body_velocity, collision_normal) - - # Coefficient of restitution (0 < e <= 1) - # Tungsten Carbide = 0.7 to 1.0 - # Steel = 0.6 to 0.9 - # High Grade aluminium hardened = 0.7 to 0.8 - # Low Grade aluminium unhardened = 0.3 to 0.4 - # Pure aluminium = 0.08 to 0.12 - restitution_coefficient = 0.5 - # println("initial_velocity_normal", initial_velocity_normal, " max ", max_overlap) - # todo: this actually depends on the timestep - if norm(initial_velocity_normal) < 1e-2 || max_overlap < 1e-6 - restitution_coefficient = 0.0 - end - - # Calculate the required change in velocity along the normal (reversing it to simulate a bounce) - velocity_change_normal = -(1 + restitution_coefficient) * initial_velocity_normal * - collision_normal - - # println("collision_normal ", collision_normal) - - post_collision_velocity = body_velocity + velocity_change_normal - - # velocity_change = norm(v_particle_system[:, 1]) - norm(v_particle_system[:, 1] + velocity_change_normal) - # if abs(velocity_change) > 1 - # println("before: ", v_particle_system[:, 1]) - # println("after: ", v_particle_system[:, 1] + velocity_change_normal) - # println("imp: ", velocity_change_normal) - - # exit(-1) - # end - - # Calculate back movement time based on post-collision velocity - back_movement_time = (max_overlap + sqrt(eps())) / norm(body_velocity) - back_movement = post_collision_velocity * back_movement_time - - # positional_correction = max_overlap * collision_normal - - # # we need to move in the opposite direction the same amount so times 2 - # positional_correction = (2.0 * max_overlap + sqrt(eps())) * collision_normal - # # if we don't move back far enough the current version gets unstable - # movement_scale = max(1.0, 0.25 * particle_radius / norm(positional_correction)) - # positional_correction_scaled = positional_correction * movement_scale - - # position_correction = (contact_point - largest_overlap_coord) - # println("position_correction", position_correction) - # println("largest_overlap_coord", largest_overlap_coord) - # println("contact_point", contact_point) - # println("post_collision_velocity", post_collision_velocity) - # println("body_velocity", body_velocity) - # exit(-1) - - # if norm(position_correction) > 0.1 - # exit(-1) - # end - - position_correction = (contact_point - largest_overlap_coord) + - back_movement - # collision_impulse = velocity_change_normal - # println("position_correction", position_correction) - # println("back_movement", back_movement) - - # Apply the change uniformly across all particles - for particle in each_moving_particle(particle_system) - v_particle_system[:, particle] = post_collision_velocity - u_particle_system[:, particle] .+= position_correction - # if min_y > u_particle_system[2, particle] - # min_y = u_particle_system[2, particle] - # end - # u_particle_system[:, particle] += [0.0, 0.25*neighbor_radius] - end - # println("miny ", min_y) +function nhs_coords(system::FluidSystem, + neighbor::BoundarySPHSystem, u) + if neighbor.ismoving[1] + return current_coordinates(u, neighbor) end - return v_particle_system -end - -# collision influences both systems -@inline function collision_interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search, - particle_system::Union{RigidSPHSystem, - TotalLagrangianSPHSystem}, - neighbor_system::Union{RigidSPHSystem, - TotalLagrangianSPHSystem}) - return dv -end - -# NHS updates -# All systems that always move update every time -function nhs_coords(system::FluidSystem, - neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, RigidSPHSystem}, - u) - return current_coordinates(u, neighbor) + # Don't update + return nothing end function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::Union{FluidSystem, RigidSPHSystem}, u) + neighbor::FluidSystem, u) return current_coordinates(u, neighbor) end -function nhs_coords(system::RigidSPHSystem, - neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, u) - return current_coordinates(u, neighbor) +function nhs_coords(system::TotalLagrangianSPHSystem, + neighbor::TotalLagrangianSPHSystem, u) + # Don't update + return nothing end -# Only update when moving -function nhs_coords(system::Union{FluidSystem, TotalLagrangianSPHSystem, RigidSPHSystem}, +function nhs_coords(system::TotalLagrangianSPHSystem, neighbor::BoundarySPHSystem, u) if neighbor.ismoving[1] return current_coordinates(u, neighbor) @@ -917,8 +620,7 @@ function nhs_coords(system::Union{FluidSystem, TotalLagrangianSPHSystem, RigidSP end function nhs_coords(system::BoundarySPHSystem, - neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, RigidSPHSystem}, u) + neighbor::FluidSystem, u) # Don't update return nothing end @@ -928,12 +630,14 @@ function nhs_coords(system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, return current_coordinates(u, neighbor) end -function nhs_coords(system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem, u) +function nhs_coords(system::BoundarySPHSystem, + neighbor::TotalLagrangianSPHSystem, u) # Don't update return nothing end -function nhs_coords(system::RigidSPHSystem, neighbor::RigidSPHSystem, u) +function nhs_coords(system::BoundarySPHSystem, + neighbor::BoundarySPHSystem, u) # Don't update return nothing end @@ -959,8 +663,7 @@ function check_configuration(boundary_system::BoundarySPHSystem, systems) end end -function check_configuration(system::Union{TotalLagrangianSPHSystem, RigidSPHSystem}, - systems) +function check_configuration(system::TotalLagrangianSPHSystem, systems) (; boundary_model) = system foreach_system(systems) do neighbor diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index dfb642eb0..76da93cf0 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -324,13 +324,34 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - nhs = get_neighborhood_search(system, neighbor_system, semi) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, v, - v_neighbor_system, nhs) + # This is an optimization for simulations with large and complex boundaries. + # Especially, in 3D simulations with large and/or complex structures outside + # of areas with permanent flow. + # Note: The version iterating neighbors first is not thread parallelizable. + # The factor is based on the achievable speed-up of the thread parallelizable version. + if nparticles(system) > + ceil(Int, 0.5 * Threads.nthreads()) * nparticles(neighbor_system) + nhs = get_neighborhood_search(neighbor_system, system, semi) + + # Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries + adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system, + system_coords, neighbor_coords, + v_neighbor_system, nhs) + else + nhs = get_neighborhood_search(system, neighbor_system, semi) + + # Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries + adami_pressure_extrapolation!(boundary_model, system, neighbor_system, + system_coords, neighbor_coords, + v_neighbor_system, nhs) + end + for particle in eachparticle(system) + # Limit pressure to be non-negative to avoid attractive forces between fluid and + # boundary particles at free surfaces (sticking artifacts). + pressure[particle] = max(pressure[particle], 0.0) + end end @trixi_timeit timer() "inverse state equation" @threaded for particle in eachparticle(system) @@ -365,9 +386,38 @@ function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZe return boundary_model end +@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, + v_neighbor_system, + neighborhood_search) + (; pressure, cache, viscosity) = boundary_model + + for_particle_neighbor(neighbor_system, system, + neighbor_coords, system_coords, + neighborhood_search; + particles=eachparticle(neighbor_system), + parallel=false) do neighbor, particle, + pos_diff, distance + # Since neighbor and particle are switched + pos_diff = -pos_diff + adami_pressure_inner!(boundary_model, system, neighbor_system::FluidSystem, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure) + end +end + +@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, + v_neighbor_system, + neighborhood_search) + return boundary_model +end + @inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system::FluidSystem, - system_coords, neighbor_coords, v_particle_system, + system_coords, neighbor_coords, v_neighbor_system, neighborhood_search) (; pressure, cache, viscosity) = boundary_model @@ -377,49 +427,40 @@ end neighborhood_search; particles=eachparticle(system)) do particle, neighbor, pos_diff, distance - density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) - kernel_weight = smoothing_kernel(boundary_model, distance) - - # if system isa RigidSPHSystem - # resulting_acc = neighbor_system.acceleration - 100.0 * system.acceleration - # else - resulting_acc = neighbor_system.acceleration - - current_acceleration(system, particle) - # end - - if system isa RigidSPHSystem - # pressure[particle] += 4000 * kernel_weight - pressure[particle] += (10000 + 0.5 * density_neighbor * norm(v_particle_system[1:ndims(system),1])^2) * kernel_weight - # println("old ", 0.5* density_neighbor *norm(v_particle_system[:,1])^2 * kernel_weight, " ", 4000 * kernel_weight) - # println("new ", 0.5* density_neighbor *norm(v_particle_system[1:ndims(system),1]- v_neighbor_system[1:ndims(system),neighbor])^2 * kernel_weight, " ", 4000 * kernel_weight) - # else - # pressure[particle] += (0.5 * density_neighbor * norm(v_neighbor_system[1:ndims(system),neighbor])^2) * kernel_weight - end - pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, - neighbor) + - dot(resulting_acc, density_neighbor * pos_diff)) * - kernel_weight - - cache.volume[particle] += kernel_weight - - # why is this done here? - compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, - kernel_weight, particle, neighbor) - end - - for particle in eachparticle(system) - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). - pressure[particle] = max(pressure[particle], 0.0) + adami_pressure_inner!(boundary_model, system, neighbor_system, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure) end end @inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, v_particle_system, + system_coords, neighbor_coords, v_neighbor_system, neighborhood_search) return boundary_model end +@inline function adami_pressure_inner!(boundary_model, system, + neighbor_system::FluidSystem, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure) + density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) + + resulting_acc = neighbor_system.acceleration - + current_acceleration(system, particle) + + kernel_weight = smoothing_kernel(boundary_model, distance) + + pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, + neighbor) + + dot(resulting_acc, density_neighbor * pos_diff)) * + kernel_weight + + cache.volume[particle] += kernel_weight + + compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, + kernel_weight, particle, neighbor) +end + function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) return cache diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 574e832d1..965e55062 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -1,8 +1,7 @@ # Interaction of boundary with other systems function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::Union{BoundarySPHSystem, RigidSPHSystem}, - neighbor_system) + particle_system::BoundarySPHSystem, neighbor_system) # TODO Solids and moving boundaries should be considered in the continuity equation return dv end @@ -10,8 +9,7 @@ end # For dummy particles with `ContinuityDensity`, solve the continuity equation function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::Union{BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, - RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}}, + particle_system::BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, neighbor_system::FluidSystem) (; boundary_model) = particle_system fluid_density_calculator = neighbor_system.density_calculator diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 1e5a325aa..b30d8fc41 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -17,11 +17,9 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM boundary_model :: BM movement :: M ismoving :: Vector{Bool} - particle_spacing :: ELTYPE cache :: C - function BoundarySPHSystem(initial_condition, model; movement=nothing, - particle_spacing=NaN) + function BoundarySPHSystem(initial_condition, model; movement=nothing) coordinates = copy(initial_condition.coordinates) NDIMS = size(coordinates, 1) ismoving = zeros(Bool, 1) @@ -37,7 +35,7 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIM return new{typeof(model), NDIMS, eltype(coordinates), typeof(movement), typeof(cache)}(initial_condition, coordinates, model, movement, - ismoving, particle_spacing, cache) + ismoving, cache) end end @@ -101,7 +99,6 @@ function Base.show(io::IO, system::BoundarySPHSystem) print(io, "BoundarySPHSystem{", ndims(system), "}(") print(io, system.boundary_model) print(io, ", ", system.movement) - print(io, ", ", system.particle_spacing) print(io, ") with ", nparticles(system), " particles") end @@ -114,7 +111,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) summary_header(io, "BoundarySPHSystem{$(ndims(system))}") summary_line(io, "#particles", nparticles(system)) summary_line(io, "boundary model", system.boundary_model) - summary_line(io, "particle spacing", system.particle_spacing) summary_line(io, "movement function", isnothing(system.movement) ? "nothing" : string(system.movement.movement_function)) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 21ed2d747..0b913bc23 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -153,40 +153,3 @@ end return p_a, p_a end - -# # Fluid particles collide with a boundary particle at distance =< 0 to the surface boundary -# # particle to prevent them going through the first layer of boundary particles -# @inline function collision_interact!(dv, v_particle_system, u_particle_system, -# v_neighbor_system, u_neighbor_system, -# neighborhood_search, particle_system::FluidSystem, -# neighbor_system::Union{BoundarySPHSystem, -# RigidSPHSystem, -# TotalLagrangianSPHSystem}) -# # (particle_spacing) = neighbor_system - -# system_coords = current_coordinates(u_particle_system, particle_system) -# neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - -# for_particle_neighbor(particle_system, neighbor_system, -# system_coords, neighbor_system_coords, -# neighborhood_search) do particle, neighbor, pos_diff, distance -# if distance <= sqrt(eps()) -# # todo: not correct -# normal = pos_diff / distance - -# part_v = extract_svector(v_particle_system, particle_system, particle) -# nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) -# # Relative velocity -# rel_vel_normal = dot(part_v - nghbr_v, normal) - -# collision_damping_coefficient = 0.1 -# force_magnitude = collision_damping_coefficient * rel_vel_normal -# force = force_magnitude * normal - -# @inbounds for i in 1:ndims(particle_system) -# dv[i, particle] += force[i] / mass[particle] -# end -# end -# end -# return dv -# end diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index f7fcc3610..5eabef80e 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -2,12 +2,7 @@ # interactions between the different system types. include("fluid/fluid.jl") include("boundary/boundary.jl") - -# needs to be after boundary.jl -include("solid/rigid_body_sph/rigid_body_sph.jl") -# needs to be after boundary.jl include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") - # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem` # and the `TotalLagrangianSPHSystem`. include("boundary/monaghan_kajtar/monaghan_kajtar.jl") @@ -17,4 +12,3 @@ include("fluid/weakly_compressible_sph/rhs.jl") include("fluid/entropically_damped_sph/rhs.jl") include("boundary/rhs.jl") include("solid/total_lagrangian_sph/rhs.jl") -include("solid/rigid_body_sph/rhs.jl") diff --git a/src/schemes/solid/rigid_body_sph/rhs.jl b/src/schemes/solid/rigid_body_sph/rhs.jl deleted file mode 100644 index acc2bd886..000000000 --- a/src/schemes/solid/rigid_body_sph/rhs.jl +++ /dev/null @@ -1,118 +0,0 @@ - -# Solid-fluid interaction -# Since the object is rigid the movement gets averaged across all particles. -function interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::RigidSPHSystem, - neighbor_system::FluidSystem) - sound_speed = system_sound_speed(neighbor_system) - - system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - total_dv = zeros(ndims(particle_system)) - - # Loop over all pairs of particles and neighbors within the kernel cutoff. - for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0. - distance < sqrt(eps()) && return - - # Apply the same force to the solid particle - # that the fluid particle experiences due to the soild particle. - # Note that the same arguments are passed here as in fluid-solid interact!, - # except that pos_diff has a flipped sign. - # - # In fluid-solid interaction, use the "hydrodynamic mass" of the solid particles - # corresponding to the rest density of the fluid and not the material density. - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 - - # Use kernel from the fluid system in order to get the same force here in - # solid-fluid interaction as for fluid-solid interaction. - grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance) - - # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles - # corresponding to the chosen boundary model. - p_a = particle_pressure(v_particle_system, particle_system, particle) - p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) - - # Particle and neighbor (and corresponding systems and all corresponding quantities) - # are switched in the following two calls. - # This way, we obtain the exact same force as for the fluid-solid interaction, - # but with a flipped sign (because `pos_diff` is flipped compared to fluid-solid). - dv_boundary = pressure_acceleration(neighbor_system, particle_system, particle, - m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, - distance, grad_kernel, - neighbor_system.correction) - - dv_viscosity_ = dv_viscosity(neighbor_system, particle_system, - v_neighbor_system, v_particle_system, - neighbor, particle, pos_diff, distance, - sound_speed, m_b, m_a, rho_mean) - - dv_particle = dv_boundary + dv_viscosity_ - - for i in 1:ndims(particle_system) - # Multiply `dv` (acceleration on fluid particle b) by the mass of - # particle b to obtain the same force as for the fluid-solid interaction. - # Divide by the material mass of particle a to obtain the acceleration - # of solid particle a. - total_dv[i] += dv[i, particle] + - dv_particle[i] * m_b / particle_system.mass[particle] - end - - continuity_equation!(dv, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) - end - - #total_dv ./=nparticles(particle_system) - for particle in each_moving_particle(particle_system) - for i in 1:ndims(particle_system) - dv[i, particle] += total_dv[i] / nparticles(particle_system) - end - end - return dv -end - -@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::RigidSPHSystem, - neighbor_system::FluidSystem, - grad_kernel) - return dv -end - -@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, - neighbor_system::FluidSystem, - grad_kernel) - fluid_density_calculator = neighbor_system.density_calculator - - v_diff = current_velocity(v_particle_system, particle_system, particle) - - current_velocity(v_neighbor_system, neighbor_system, neighbor) - - # Call the dummy BC version of the continuity equation - continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff, - grad_kernel, particle) -end - -# Solid-boundary interaction -function interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::RigidSPHSystem, - neighbor_system::Union{BoundarySPHSystem}) - - # TODO continuity equation? - return dv -end diff --git a/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl b/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl deleted file mode 100644 index c67930ab8..000000000 --- a/src/schemes/solid/rigid_body_sph/rigid_body_sph.jl +++ /dev/null @@ -1 +0,0 @@ -include("system.jl") diff --git a/src/schemes/solid/rigid_body_sph/system.jl b/src/schemes/solid/rigid_body_sph/system.jl deleted file mode 100644 index 44de55356..000000000 --- a/src/schemes/solid/rigid_body_sph/system.jl +++ /dev/null @@ -1,205 +0,0 @@ -@doc raw""" - RigidSPHSystem(initial_condition; - boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, NDIMS)) - -System for a SPH particles of a rigid structure. - -# Arguments -- `initial_condition`: Initial condition representing the system's particles. - -# Keyword Arguments -- `boundary_model`: Boundary model to compute the hydrodynamic density and pressure for - fluid-structure interaction (see [Boundary Models](@ref boundary_models)). -- `acceleration`: Acceleration vector for the system. (default: zero vector) -""" -struct RigidSPHSystem{BM, NDIMS, ELTYPE <: Real} <: SolidSystem{NDIMS} - initial_condition :: InitialCondition{ELTYPE} - local_coordinates :: Array{ELTYPE, 2} # [dimension, particle] - mass :: Array{ELTYPE, 1} # [particle] - material_density :: Array{ELTYPE, 1} # [particle] - acceleration :: SVector{NDIMS, ELTYPE} - center_of_gravity :: Array{ELTYPE, 1} # [dimension] - collision_impulse :: Array{ELTYPE, 1} # [dimension] - collision_u :: Array{ELTYPE, 1} # [dimension] - particle_spacing :: ELTYPE - has_collided :: MutableBool - boundary_model :: BM - - function RigidSPHSystem(initial_condition; boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, ndims(initial_condition)), - particle_spacing=NaN) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - # Make acceleration an SVector - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end - - local_coordinates = copy(initial_condition.coordinates) - mass = copy(initial_condition.mass) - material_density = copy(initial_condition.density) - - # TODO: calculate center of gravity - cog = zeros(SVector{NDIMS, ELTYPE}) - collision_impulse = zeros(SVector{NDIMS, ELTYPE}) - collision_u = zeros(SVector{NDIMS, ELTYPE}) - - return new{typeof(boundary_model), NDIMS, ELTYPE}(initial_condition, - local_coordinates, - mass, material_density, - acceleration_, cog, - collision_impulse, - collision_u, - particle_spacing, - MutableBool(false), - boundary_model) - end -end - -# Info output functions - -function Base.show(io::IO, system::RigidSPHSystem) - @nospecialize system # reduce precompilation time - - print(io, "RigidSPHSystem{", ndims(system), "}(") - print(io, ", ", system.acceleration) - print(io, ", ", system.boundary_model) - print(io, ") with ", nparticles(system), " particles") -end - -function Base.show(io::IO, ::MIME"text/plain", system::RigidSPHSystem) - @nospecialize system # reduce precompilation time - - if get(io, :compact, false) - show(io, system) - else - summary_header(io, "RigidSPHSystem{$(ndims(system))}") - summary_line(io, "total #particles", nparticles(system)) - summary_line(io, "acceleration", system.acceleration) - summary_line(io, "boundary model", system.boundary_model) - summary_footer(io) - end -end - -# Memory allocation accessors - -@inline function v_nvariables(system::RigidSPHSystem) - return ndims(system) -end - -@inline function v_nvariables(system::RigidSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}) - return ndims(system) + 1 -end - -# Variable accessors - -@inline function n_moving_particles(system::RigidSPHSystem) - return nparticles(system) -end - -@inline local_coordinates(system::RigidSPHSystem) = system.local_coordinates - -# @inline function current_coordinates(u, system::RigidSPHSystem) -# return system.coordinates -# end - -@inline function current_velocity(v, system::RigidSPHSystem, particle) - return extract_svector(v, system, particle) -end - -@inline function viscous_velocity(v, system::RigidSPHSystem, particle) - return extract_svector(system.boundary_model.cache.wall_velocity, system, particle) -end - -@inline function particle_density(v, system::RigidSPHSystem, particle) - return particle_density(v, system.boundary_model, system, particle) -end - -# In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles -# corresponding to the chosen boundary model. -@inline function particle_pressure(v, system::RigidSPHSystem, particle) - return particle_pressure(v, system.boundary_model, system, particle) -end - -@inline function hydrodynamic_mass(system::RigidSPHSystem, particle) - return system.boundary_model.hydrodynamic_mass[particle] -end - -function initialize!(system::RigidSPHSystem, neighborhood_search) -end - -function write_u0!(u0, system::RigidSPHSystem) - (; initial_condition) = system - - for particle in each_moving_particle(system) - # Write particle coordinates - for dim in 1:ndims(system) - u0[dim, particle] = initial_condition.coordinates[dim, particle] - end - end - - return u0 -end - -function write_v0!(v0, system::RigidSPHSystem) - (; initial_condition, boundary_model) = system - - for particle in each_moving_particle(system) - # Write particle velocities - for dim in 1:ndims(system) - v0[dim, particle] = initial_condition.velocity[dim, particle] - end - end - - write_v0!(v0, boundary_model, system) - - return v0 -end - -function write_v0!(v0, model, system::RigidSPHSystem) - return v0 -end - -function write_v0!(v0, ::BoundaryModelDummyParticles{ContinuityDensity}, - system::RigidSPHSystem) - (; cache) = system.boundary_model - (; initial_density) = cache - - for particle in each_moving_particle(system) - # Set particle densities - v0[ndims(system) + 1, particle] = initial_density[particle] - end - - return v0 -end - -# Update functions -# function update_positions!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) -# # (; movement) = system - -# # movement(system, t) -# end - -# function update_quantities!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) -# return system -# end - -function update_final!(system::RigidSPHSystem, v, u, v_ode, u_ode, semi, t) - (; boundary_model) = system - - # Only update boundary model - update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) -end - -function viscosity_model(system::RigidSPHSystem) - return system.boundary_model.viscosity -end - -function calculate_dt(v_ode, u_ode, cfl_number, system::RigidSPHSystem) - (; acceleration, particle_spacing) = system - - return cfl_number * particle_spacing/maximum(abs.(acceleration)) -end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index b7e7f3f6f..9db300269 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -252,13 +252,6 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end -function write2vtk!(vtk, v, u, t, system::RigidSPHSystem; write_meta_data=true) - vtk["velocity"] = view(v, 1:ndims(system), :) - vtk["material_density"] = system.material_density - - write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) -end - function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true) write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl index 3076b9274..0a97850bd 100644 --- a/test/callbacks/callbacks.jl +++ b/test/callbacks/callbacks.jl @@ -2,4 +2,5 @@ include("info.jl") include("stepsize.jl") include("postprocess.jl") + include("solution_saving.jl") end diff --git a/test/callbacks/solution_saving.jl b/test/callbacks/solution_saving.jl new file mode 100644 index 000000000..3cc765b65 --- /dev/null +++ b/test/callbacks/solution_saving.jl @@ -0,0 +1,68 @@ +@testset verbose=true "SolutionSavingCallback" begin + @testset verbose=true "show" begin + out = joinpath(pkgdir(TrixiParticles), "out") + output_directory_padded = out * " "^(65 - length(out)) + + @testset verbose=true "dt" begin + callback = SolutionSavingCallback(dt=0.02, prefix="test", output_directory=out) + + show_compact = "SolutionSavingCallback(dt=0.02)" + @test repr(callback) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ SolutionSavingCallback │ + │ ══════════════════════ │ + │ dt: ……………………………………………………………………… 0.02 │ + │ custom quantities: ……………………………… nothing │ + │ save initial solution: …………………… yes │ + │ save final solution: ………………………… yes │ + │ output directory: ………………………………… $(output_directory_padded)│ + │ prefix: …………………………………………………………… test │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + + @testset verbose=true "interval" begin + callback = SolutionSavingCallback(interval=100, prefix="test", + output_directory=out) + + show_compact = "SolutionSavingCallback(interval=100)" + @test repr(callback) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ SolutionSavingCallback │ + │ ══════════════════════ │ + │ interval: ……………………………………………………… 100 │ + │ custom quantities: ……………………………… nothing │ + │ save initial solution: …………………… yes │ + │ save final solution: ………………………… yes │ + │ output directory: ………………………………… $(output_directory_padded)│ + │ prefix: …………………………………………………………… test │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + + @testset verbose=true "interval" begin + callback = SolutionSavingCallback(save_times=[1.0, 2.0, 3.0], prefix="test", + output_directory=out) + + show_compact = "SolutionSavingCallback(save_times=[1.0, 2.0, 3.0])" + @test repr(callback) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ SolutionSavingCallback │ + │ ══════════════════════ │ + │ save_times: ………………………………………………… [1.0, 2.0, 3.0] │ + │ custom quantities: ……………………………… nothing │ + │ save initial solution: …………………… yes │ + │ save final solution: ………………………… yes │ + │ output directory: ………………………………… $(output_directory_padded)│ + │ prefix: …………………………………………………………… test │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + end +end diff --git a/test/validation/validation.jl b/test/validation/validation.jl index 95dcfe1bd..97413c710 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -46,16 +46,18 @@ ] @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 - @test isapprox(error_edac_P1, 0, atol=1e-9) - @test isapprox(error_edac_P2, 0, atol=1e-11) if VERSION >= v"1.10" + @test isapprox(error_edac_P1, 0, atol=eps()) + @test isapprox(error_edac_P2, 0, atol=eps()) @test isapprox(error_wcsph_P1, 0, atol=eps()) @test isapprox(error_wcsph_P2, 0, atol=eps()) else # 1.9 causes a large difference in the solution - @test isapprox(error_wcsph_P1, 0, atol=10) - @test isapprox(error_wcsph_P2, 0, atol=1e-3) + @test isapprox(error_edac_P1, 0, atol=eps()) + @test isapprox(error_edac_P2, 0, atol=eps()) + @test isapprox(error_wcsph_P1, 0, atol=eps()) + @test isapprox(error_wcsph_P2, 0, atol=eps()) end # Ignore method redefinitions from duplicate `include("../validation_util.jl")` diff --git a/validation/dam_break_2d/validation_dam_break_2d.jl b/validation/dam_break_2d/validation_dam_break_2d.jl index bfe216fde..4bf57d59e 100644 --- a/validation/dam_break_2d/validation_dam_break_2d.jl +++ b/validation/dam_break_2d/validation_dam_break_2d.jl @@ -115,12 +115,12 @@ run_file_edac_name = joinpath("out", reference_data = JSON.parsefile(reference_file_edac_name) run_data = JSON.parsefile(run_file_edac_name) -error_edac_P1 = interpolated_mse(reference_data["pressure_P1_fluid_1"]["time"], +error_edac_P1 = interpolated_mre(reference_data["pressure_P1_fluid_1"]["time"], reference_data["pressure_P1_fluid_1"]["values"], run_data["pressure_P1_fluid_1"]["time"], run_data["pressure_P1_fluid_1"]["values"]) -error_edac_P2 = interpolated_mse(reference_data["pressure_P2_fluid_1"]["time"], +error_edac_P2 = interpolated_mre(reference_data["pressure_P2_fluid_1"]["time"], reference_data["pressure_P2_fluid_1"]["values"], run_data["pressure_P2_fluid_1"]["time"], run_data["pressure_P2_fluid_1"]["values"]) @@ -150,12 +150,12 @@ run_file_wcsph_name = joinpath("out", reference_data = JSON.parsefile(reference_file_wcsph_name) run_data = JSON.parsefile(run_file_wcsph_name) -error_wcsph_P1 = interpolated_mse(reference_data["pressure_P1_fluid_1"]["time"], +error_wcsph_P1 = interpolated_mre(reference_data["pressure_P1_fluid_1"]["time"], reference_data["pressure_P1_fluid_1"]["values"], run_data["pressure_P1_fluid_1"]["time"], run_data["pressure_P1_fluid_1"]["values"]) -error_wcsph_P2 = interpolated_mse(reference_data["pressure_P2_fluid_1"]["time"], +error_wcsph_P2 = interpolated_mre(reference_data["pressure_P2_fluid_1"]["time"], reference_data["pressure_P2_fluid_1"]["values"], run_data["pressure_P2_fluid_1"]["time"], run_data["pressure_P2_fluid_1"]["values"]) diff --git a/validation/dam_break_2d/validation_reference_edac_0015.json b/validation/dam_break_2d/validation_reference_edac_0015.json index 24757b07f..fb8ead957 100644 --- a/validation/dam_break_2d/validation_reference_edac_0015.json +++ b/validation/dam_break_2d/validation_reference_edac_0015.json @@ -139,79 +139,79 @@ 0.0, 0.0, 0.0, - 604.3957582215899, - 854.3064249696633, - 806.7433554643078, - 863.2834525307577, - 1172.7188760376416, - 1185.7354823999622, - 1029.8194599143999, - 1220.8424591423286, - 1129.1130583244467, - 1681.2917027879928, - 1541.8470498978272, - 1941.6817449062, - 1703.3428527373653, - 1820.7752456908524, - 2273.61907013479, - 2059.549527721113, - 2492.4265793163177, - 2431.891436605787, - 2341.4127018992012, - 2380.505662332444, - 2476.6304162505776, - 2377.2152352341263, - 2407.6934010063724, - 2621.6411341900084, - 2604.6982173831625, - 2546.399445810557, - 2689.684321946729, - 2654.6394300204142, - 2815.3615097824277, - 2642.321188486082, - 2862.2882426968963, - 2811.52377372663, - 2800.9092542116946, - 2952.956412783148, - 2844.263845334785, - 2874.303528967565, - 2947.9713667109695, - 2981.7363562084547, - 3022.715296426873, - 3095.8959184057803, - 3223.2834166674697, - 3438.1549404217467, - 3461.6157573942837, - 3889.5328111053327, - 4764.835081898881, - 6691.330625903637, - 5937.573914006671, - 6470.044187005868, - 5114.121150798657, - 3634.979008335772, - 3802.039385890102, - 2751.9598060963713, - 2027.890808792743, - 1718.7065304962896, - 1535.5876413862059, - 1382.3526450110292, - 1479.4637909212727, - 1636.5432932362105, - 1580.8409259586185, - 1362.2791052129037, - 1307.4981915271765, - 1483.094125847357, - 1582.6083889647803, - 3052.2515396988983, - 34536.296819296316, - 37.65056335694383 + 604.3957582210057, + 854.3064249701162, + 806.7433554645238, + 863.2834525323564, + 1172.7188760362817, + 1185.7354823957364, + 1029.8194599051196, + 1220.8424591289627, + 1129.1130583250354, + 1681.2917027821925, + 1541.84704989596, + 1941.681744886013, + 1703.342852837277, + 1820.7752456709127, + 2273.6190702267386, + 2059.549528549632, + 2492.42657872154, + 2431.891435051085, + 2341.412705521953, + 2380.5056651594596, + 2476.6304064185897, + 2377.215248491543, + 2407.6934046918395, + 2621.6411290967458, + 2604.698214284965, + 2546.399443630783, + 2689.6843234734833, + 2654.639435738155, + 2815.3615123943055, + 2642.321187463903, + 2862.2882422981534, + 2811.523774221353, + 2800.9092531867764, + 2952.9564218334053, + 2844.26385035426, + 2874.3035321104517, + 2947.9713669965477, + 2981.736356294694, + 3022.7152925706846, + 3095.8959107824994, + 3223.283421806281, + 3438.1549422956728, + 3461.615757497162, + 3889.532826771071, + 4764.8350876461145, + 6691.330618320588, + 5937.57389697787, + 6470.044184482702, + 5114.121178496523, + 3634.978964640017, + 3802.0393876762696, + 2751.9597721194446, + 2027.8908055573243, + 1718.7065075201758, + 1535.5875986343772, + 1382.352588739524, + 1479.4637838388012, + 1636.543321007025, + 1580.840904760441, + 1362.2790951007676, + 1307.4981559335906, + 1483.0941976925808, + 1582.6084515968898, + 3052.2514672201683, + 34536.29685017654, + 37.65095815197997 ], "datatype": "Float64", "type": "series" }, "meta": { "julia_version": "1.10.2", - "solver_version": "dcd1d2b1-dirty", + "solver_version": "v0.1.0-22-g4d428a23-dirty", "solver_name": "TrixiParticles.jl" }, "pressure_P2_fluid_1": { @@ -360,43 +360,43 @@ 0.0, 0.0, 0.0, - 1.4586145070841525, - 35.23970220674348, - 37.54572391444544, - 13.731385167742658, - 18.846457975737682, - 19.8224841977086, - 16.230770674916435, - 56.876420772446885, - 37.33671454315734, - 78.42964736040615, - 75.04969119915044, - 166.61484088498904, - 162.7352955160456, - 281.7282268014637, - 314.52179459013166, - 297.23412867557687, - 357.6400910860976, - 448.6558207559398, - 465.072910798512, - 612.2401011562348, - 751.9411583776679, - 889.9451285845244, - 1033.4736483837655, - 1058.202206530216, - 1247.2730240894725, - 1414.0754260543076, - 1588.386535262655, - 1665.7826025761042, - 1545.8294662391331, - 1673.9628009744847, - 1683.2228289778777, - 1629.4962992741189, - 1359.6993719827465, - 887.91187263999, - 1629.9983580422795, - 1234.4717693503253, - 2004.3258881084407, + 1.4586145071541898, + 35.23970220674639, + 37.545723914493706, + 13.731385167730085, + 18.84645797568068, + 19.822484197824288, + 16.23077067519315, + 56.87642077158365, + 37.33671454346889, + 78.42964734337639, + 75.04969120367215, + 166.6148408881951, + 162.73529547483014, + 281.72822678093445, + 314.5217962647108, + 297.23412415145486, + 357.6400963258884, + 448.65581898807187, + 465.0729072685026, + 612.2400999708409, + 751.9411576661471, + 889.945132047711, + 1033.4736519596122, + 1058.2021923062105, + 1247.2730353264728, + 1414.0754183693182, + 1588.3865333590663, + 1665.7825902521158, + 1545.8294821241884, + 1673.9628180445059, + 1683.2228111190875, + 1629.4962936654842, + 1359.6993745648126, + 887.9118718296219, + 1629.9983551561381, + 1234.4717618954708, + 2004.325882177016, 0.0, 0.0, 0.0, @@ -538,98 +538,98 @@ 1.302960573952499, 1.3475843440024642, 1.3976361698571822, - 1.4522199811561385, - 1.510520947475266, - 1.57214293886014, - 1.6367982265986314, - 1.7042518338466675, - 1.7742937998709938, - 1.8467352027968718, - 1.9213758913869927, - 1.9980351614993204, - 2.076544579209514, - 2.1566383641391065, - 2.238003134643121, - 2.3206381687170157, - 2.4044919394042736, - 2.489540630902278, - 2.5757273488305, - 2.6627522454913684, - 2.7501950551381626, - 2.837835825711921, - 2.9257131544255777, - 3.013897557214555, - 3.102445768668488, - 3.1914847216832074, - 3.2105882666919, - 3.2106025664774633, - 3.210407075423367, - 3.210274207277424, - 3.210140068534914, - 3.210439644792984, - 3.2112884336091394, - 3.2121208254868923, - 3.2125141820856236, - 3.212599029590728, - 3.2119451344617342, - 3.210893702477008, - 3.2095320951690165, - 3.208866042695388, - 3.20871345061095, - 3.208573760872467, - 3.2083568996870384, - 3.208389043560661, - 3.208031877949919, - 3.2079261547153943, - 3.207855507755997, - 3.207774119588862, - 3.207706002731428, - 3.2076706691657413, + 1.452219981156139, + 1.5105209474752668, + 1.572142938860141, + 1.6367982265986325, + 1.7042518338466686, + 1.774293799870995, + 1.8467352027968729, + 1.9213758913869938, + 1.9980351614993215, + 2.0765445792095147, + 2.1566383641391074, + 2.238003134643122, + 2.3206381687170166, + 2.404491939404275, + 2.4895406309022796, + 2.575727348830502, + 2.662752245491369, + 2.7501950551381618, + 2.837835825711919, + 2.925713154425575, + 3.0138975572145523, + 3.102445768668486, + 3.191484721683206, + 3.2105882666919006, + 3.210602566477464, + 3.2104070754233653, + 3.2102742072774215, + 3.210140068534911, + 3.210439644793008, + 3.211288433609166, + 3.21212082548689, + 3.2125141820856595, + 3.212599029590776, + 3.211945134461789, + 3.210893702477077, + 3.2095320951690947, + 3.2088660426953854, + 3.2087134506109494, + 3.20857376087247, + 3.2083568996870397, + 3.2083890435617195, + 3.2080318779499213, + 3.2079261547153983, + 3.207855507756003, + 3.2077741195888656, + 3.2077060027314293, + 3.207670669165742, 3.2076288978754195, - 3.2075941374177184, + 3.2075941374177193, 3.2075586504257148, - 3.2075368319522997, - 3.2075265213627904, - 3.210218793361137, - 3.21269507520042, - 3.2119664146378195, - 3.2074630303747567, - 3.2074452874588837, - 3.2074252902484437, - 3.2074042323916974, - 3.207389539186995, - 3.207376753579591, - 3.2073732467888747, - 3.207348313206412, - 3.2076648180581375, - 3.2073154618525073, - 3.2078540055088123, - 3.2085800521157415, - 3.2072818774036778, - 3.2072679757668348, - 3.207248958752789, - 3.2072243341066833, - 3.2072195892028104, - 3.2072071351452363, - 3.207191617668054, - 3.207167720999114, - 3.207148223076407, - 3.207125308504274, - 3.2071076627463984, - 3.2070796885516786, - 3.2070631437453265, - 3.207023980424339, - 3.207002203054955, - 3.2069578815542887, - 3.206882240708725, - 3.206818992733726, - 3.2067228411656616, - 3.2066412586740536, - 3.206517043308112, - 3.2063950957148974, - 3.206329455325393, - 3.206229289248847, - 3.2061593054711017 + 3.207536831952301, + 3.2075265213627917, + 3.2102187933610122, + 3.2126950751999495, + 3.21196641463669, + 3.2074630303747598, + 3.207445287458887, + 3.2074252902484464, + 3.2074042323917, + 3.2073895391869978, + 3.207376753579594, + 3.2073732467888774, + 3.2073483132064147, + 3.207664818066874, + 3.20731546185251, + 3.207854005492601, + 3.208580052098399, + 3.2072818774036804, + 3.207267975766838, + 3.2072489587527926, + 3.2072243341066873, + 3.207219589202814, + 3.207207135145238, + 3.207191617668057, + 3.207167720999116, + 3.2071482230764086, + 3.2071253085042755, + 3.2071076627464006, + 3.207079688551678, + 3.2070631437453248, + 3.207023980424338, + 3.2070022030549534, + 3.2069578815542865, + 3.206882240708722, + 3.206818992733722, + 3.206722841165657, + 3.206641258674049, + 3.20651704330811, + 3.2063950957148926, + 3.206329455325389, + 3.206229289248844, + 3.2061593054711004 ], "datatype": "Float64", "type": "series" @@ -792,16 +792,16 @@ 0.0, 0.0, 0.0, - 3.9648485713304686, - 5.688115787444037, - 5.021169068151417, - 3.4690420734565928, - 3.558737480081865, - 3.1571728411890927, - 3.002903034691488, - 4.056831531943184, - 6.550226833801881, - 4.5706557048191865, + 3.9648485713534356, + 5.688115787337868, + 5.021169068061549, + 3.4690420744982595, + 3.5587374771061193, + 3.1571728370072245, + 3.002903045034488, + 4.056831537710453, + 6.550226835881946, + 4.570655703611319, 0.0, 0.0, 0.0, @@ -809,7 +809,7 @@ 0.0, 0.0, 0.0, - 13.988872277177816, + 13.98887227719292, 0.0, 0.0, 0.0, @@ -844,4 +844,4 @@ "datatype": "Float64", "type": "series" } -} \ No newline at end of file +} diff --git a/validation/dam_break_2d/validation_reference_wcsph_0015.json b/validation/dam_break_2d/validation_reference_wcsph_0015.json index 3a3d2221e..3b0f95b9f 100644 --- a/validation/dam_break_2d/validation_reference_wcsph_0015.json +++ b/validation/dam_break_2d/validation_reference_wcsph_0015.json @@ -139,79 +139,79 @@ 0.0, 0.0, 0.0, - 741.9939338248116, - 1074.8234704502822, - 780.1667963911615, - 750.2648649048765, - 1867.0876734094586, - 2289.2233790766354, - 1854.778595920126, - 1898.4536968839027, - 1467.7442216200166, - 1634.062582456003, - 1051.0175372038652, - 1792.992475880798, - 2224.4495761272037, - 2072.49744240541, - 1759.7942677800675, - 2213.446522638474, - 2195.8824872685705, - 2296.936590368152, - 2110.2619949482696, - 2199.2839820790064, - 2970.3733860698267, - 2306.9379237499984, - 2708.549145872953, - 3000.178544197365, - 2269.792696086254, - 2084.4626138591966, - 2584.8924415020083, - 2922.9580586157217, - 2700.961322015944, - 2912.960202391661, - 2333.093486494811, - 2272.6997429718494, - 2482.9264961596828, - 2818.1473417252982, - 2995.432486972956, - 2909.3768298891655, - 3100.77021995143, - 3025.9649439094655, - 3235.1990309725757, - 3257.7483119248636, - 2988.7945918465675, - 3115.2094695758724, - 3505.9426755013787, - 3713.8656343932016, - 4586.51617966487, - 4512.604003095575, - 6079.787376393051, - 6651.042779454037, - 5998.169882646548, - 3485.069906443089, - 810.5465351054376, - 1498.1154610674191, - 1280.9451536023262, - 6123.104036824773, - 1207.454673045513, - 1872.428135603674, - 345.5216355544876, - 1190.3488951150066, - 1160.696857170737, - 1552.7750026162616, - 2007.8328589553835, - 1229.3505526363801, - 941.7174210262225, - 1420.6537905276475, - 4080.222555513338, - 774.1134490223096 + 741.993933814252, + 1074.8234704492293, + 780.1667963845828, + 750.264864820256, + 1867.0876722351645, + 2289.2233808791984, + 1854.7785895125223, + 1898.453697957532, + 1467.7442227383228, + 1634.0625899997626, + 1051.017537937822, + 1792.9924495715281, + 2224.4495834290674, + 2072.497434035885, + 1759.7943013144645, + 2213.4464461574053, + 2195.8826026854586, + 2296.9363606455972, + 2110.2618680074606, + 2199.283897208698, + 2970.3734738765897, + 2306.9378427287784, + 2708.5493847934304, + 3000.1784782069417, + 2269.7929805781623, + 2084.4623649717955, + 2584.892142305768, + 2922.95823696328, + 2700.9623996934047, + 2912.9611466447036, + 2333.0917554209454, + 2272.7000065438688, + 2482.9258774761342, + 2818.1463732338643, + 2995.433781965522, + 2909.3736079646987, + 3100.7712408135812, + 3025.9641406700944, + 3235.2029348584956, + 3257.7410461563277, + 2988.8181504113973, + 3115.194724294498, + 3505.92247124294, + 3713.8459250716937, + 4586.526088795744, + 4512.65880424846, + 6079.986577464039, + 6651.130159524973, + 5997.838678899889, + 3484.6443631044467, + 810.4039689089993, + 1497.7738917876206, + 1281.226457082536, + 6127.035664299017, + 1199.4753269710066, + 1869.222563929746, + 345.0423599469749, + 1192.8963322296063, + 1164.7884692609991, + 1552.9964308420076, + 2007.227306493361, + 1229.1139656462394, + 940.9388974307137, + 1420.324488576742, + 4081.141427086134, + 778.984800800608 ], "datatype": "Float64", "type": "series" }, "meta": { "julia_version": "1.10.2", - "solver_version": "dcd1d2b1", + "solver_version": "v0.1.0-22-g4d428a23-dirty", "solver_name": "TrixiParticles.jl" }, "pressure_P2_fluid_1": { @@ -360,41 +360,41 @@ 0.0, 0.0, 0.0, - 47.78581694313901, - 44.38574442306553, - 30.473857197402413, - 8.477074724627595, - 33.57951875059658, - 1.8632501617021102, - 22.444845391614827, - 5.229751210510513, - 12.452368822405496, - 113.28031750326882, - 16.872533929605122, - 183.53117184170725, - 124.91979697908702, - 119.72035918066965, - 157.18678733719997, - 274.41224922147677, - 366.928403397206, - 436.7982051917158, - 656.8646154705689, - 693.6804508905344, - 709.8624537628132, - 940.2761580309083, - 713.5961440910972, - 1152.7121152728616, - 1175.0856672128707, - 1117.0740820240865, - 1473.3118337985009, - 1398.3788581542326, - 1426.7186616370407, - 1724.6680215045349, - 1689.582829806212, - 1716.4125370027475, - 2313.2602797513496, - 1548.2978703765393, - 872.8096099036991, + 47.78581701120133, + 44.38574448700575, + 30.47385719369803, + 8.477074777805942, + 33.57951839840989, + 1.8632502205569363, + 22.444846119437738, + 5.229750916569136, + 12.452367852533479, + 113.28032802403966, + 16.872523440014806, + 183.53111225036926, + 124.91979335190122, + 119.72020754524598, + 157.186868843168, + 274.4122531283022, + 366.92834558395305, + 436.79833075915656, + 656.8648337902086, + 693.6806185748894, + 709.8624031802605, + 940.2755903496779, + 713.5958561086439, + 1152.7132323277215, + 1175.0858756074435, + 1117.0732832554718, + 1473.3114994251675, + 1398.3800311373202, + 1426.7187322271566, + 1724.6666661737895, + 1689.5894191979883, + 1716.4136340314021, + 2313.2505906487645, + 1548.2452467272121, + 872.8944564735984, 0.0, 0.0, 0.0, @@ -531,105 +531,105 @@ "system_name": "fluid", "values": [ 1.1925, - 1.1979096218668643, - 1.2114046382489774, - 1.233642069008528, - 1.2644066208331137, - 1.3028362471750408, - 1.347675181367617, - 1.3979265161874306, - 1.4526822092129996, - 1.5111470319630962, - 1.5729272010563125, - 1.6376619970447728, - 1.7050782695517812, - 1.774922235352398, - 1.8470381017706874, - 1.9214486425691797, - 1.9977114877634126, - 2.0756036806364517, + 1.1979096218668641, + 1.2114046382489772, + 1.2336420690085275, + 1.264406620833113, + 1.30283624717504, + 1.3476751813676164, + 1.3979265161874295, + 1.4526822092129983, + 1.5111470319630946, + 1.572927201056311, + 1.6376619970447712, + 1.7050782695517803, + 1.7749222353523972, + 1.8470381017706896, + 1.9214486425691817, + 1.9977114877634146, + 2.075603680636452, 2.1548845384780377, - 2.235402446374903, - 2.317146330306623, - 2.4000829664557464, - 2.484138084380929, - 2.5693227139746573, - 2.6555845619869327, - 2.7426282187297657, - 2.830107717849038, - 2.9177936457239837, + 2.2354024463749025, + 2.317146330306622, + 2.400082966455745, + 2.4841380843809278, + 2.569322713974655, + 2.6555845619869305, + 2.742628218729764, + 2.8301077178490366, + 2.917793645723983, 3.0055781193133417, 3.0934656910674185, - 3.1819553106957805, - 3.2108247664460037, + 3.1819553106957756, + 3.210824766445999, 3.2113869376507562, 3.2114127786783255, - 3.2114161066075466, - 3.2113893826012587, - 3.211426416730024, - 3.2113609775292216, - 3.21134242815772, - 3.2113732404825783, - 3.2114163514270784, - 3.2115424540866804, - 3.2114240008532415, - 3.21133982171087, - 3.211435764568252, - 3.2114307576517303, - 3.2113884470340004, - 3.211430796900016, - 3.21138905570089, - 3.211385392690142, - 3.211391389367187, - 3.2113867220811616, - 3.211314609628206, - 3.2112229705963546, - 3.211144904724678, - 3.2110067641757443, - 3.210911971600086, - 3.2107740015910333, - 3.2105999166679235, - 3.210409429662026, - 3.2101913618122015, - 3.209950737465752, - 3.2096723630110047, - 3.209388535442989, - 3.209054807974806, - 3.2087189035854937, - 3.2084002150370896, - 3.2092529975678827, - 3.210961614417827, - 3.2087142145674608, - 3.2078771017362766, - 3.2099811740892332, - 3.2099475931241868, - 3.2057948825772176, - 3.2053841685407347, - 3.206091263909073, - 3.205519339413853, - 3.2060202862612424, - 3.208319007154801, - 3.2093695706812637, - 3.2093274365942297, - 3.2077697929876843, - 3.2052776717641263, - 3.204666763480492, - 3.204935394911624, - 3.205533549952792, - 3.206955192729286, - 3.205966890849804, - 3.205813036921348, - 3.2054758271429558, - 3.2056931166230225, - 3.2057049469317778, - 3.205688541125392, - 3.2056979581817227, - 3.205641713595731, - 3.207203044444701, - 3.2078958675858784, - 3.2077053880890443, - 3.2064670429651616, - 3.2020423230291524 + 3.211416106607547, + 3.2113893826012614, + 3.211426416730049, + 3.2113609775292393, + 3.2113424281577405, + 3.2113732404827138, + 3.2114163514266383, + 3.2115424540877786, + 3.2114240008557804, + 3.2113398217104034, + 3.2114357645677893, + 3.2114307576516063, + 3.2113884470339853, + 3.211430796900003, + 3.211389055700846, + 3.2113853926900813, + 3.21139138936708, + 3.211386722081144, + 3.2113146096279768, + 3.2112229705960287, + 3.2111449047243386, + 3.2110067641753868, + 3.2109119715996353, + 3.210774001590455, + 3.210599916667371, + 3.210409429661165, + 3.2101913618115168, + 3.20995073746505, + 3.209672363010088, + 3.209388535442117, + 3.2090548079737142, + 3.2087189035843537, + 3.2084002150359723, + 3.209253007341981, + 3.210961634193927, + 3.208714219483822, + 3.2078770543055657, + 3.2099811750510763, + 3.2099476105422085, + 3.205794882576252, + 3.205384168539671, + 3.20609132730366, + 3.2055193740331265, + 3.206017846407191, + 3.208320022840661, + 3.209364084855916, + 3.2093323815411345, + 3.2077746632044515, + 3.2052718199939743, + 3.204666769872325, + 3.2049353274677914, + 3.205556387328138, + 3.2069851034663706, + 3.2059668637455325, + 3.205811839035014, + 3.205474317245427, + 3.205693324445893, + 3.205705272070128, + 3.2056899532532865, + 3.205701459845998, + 3.2056470262047125, + 3.2071535188963542, + 3.207886027964543, + 3.207631231852247, + 3.2062764850292496, + 3.202043659162054 ], "datatype": "Float64", "type": "series" @@ -789,29 +789,29 @@ 0.0, 0.0, 0.0, - 7.496013435689593, - 4.405923582489511, - 3.863763140478915, - 5.608197870347884, - 5.873367887336063, - 8.085513458105599, - 5.478876868208552, - 4.599394355614402, - 4.4993063713207775, - 4.593781791458133, - 4.183019826563516, - 1.4104621904133992, - 3.2046660479543077, - 11.008713759804175, - 5.743853671430627, - 4.599785421280336, + 7.496013434489482, + 4.4059235799132725, + 3.863763129703675, + 5.608197865483601, + 5.873368216011654, + 8.085513300597075, + 5.478876262291463, + 4.599394758694961, + 4.499307194633306, + 4.593780734393691, + 4.18301483779327, + 1.4104607407747796, + 3.2046657016293922, + 11.008712667545892, + 5.743854966354466, + 4.599781912734178, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 6.279733295239109, + 6.279727857215189, 0.0, 0.0, 0.0, @@ -844,4 +844,4 @@ "datatype": "Float64", "type": "series" } -} \ No newline at end of file +} diff --git a/validation/validation_util.jl b/validation/validation_util.jl index 01db69635..56433f0f0 100644 --- a/validation/validation_util.jl +++ b/validation/validation_util.jl @@ -1,3 +1,9 @@ +# Perform linear interpolation to find a value at `interpolation_point` using arrays `x` and `y`. +# +# Arguments: +# - `x` : The array of abscissas (e.g., time points). +# - `y` : The array of ordinates (e.g., data points corresponding to `x`). +# - `interpolation_point` : The point at which to interpolate the data. function linear_interpolation(x, y, interpolation_point) if !(first(x) <= interpolation_point <= last(x)) throw(ArgumentError("`interpolation_point` at $interpolation_point is outside the interpolation range")) @@ -12,6 +18,15 @@ function linear_interpolation(x, y, interpolation_point) return y[i] + slope * (interpolation_point - x[i]) end + +# Calculate the mean squared error (MSE) between interpolated simulation values and reference values +# over a common time range. +# +# Arguments: +# - `reference_time` : Time points for the reference data. +# - `reference_values` : Data points for the reference data. +# - `simulation_time` : Time points for the simulation data. +# - `simulation_values` : Data points for the simulation data. function interpolated_mse(reference_time, reference_values, simulation_time, simulation_values) if last(simulation_time) > last(reference_time) @@ -34,6 +49,37 @@ function interpolated_mse(reference_time, reference_values, simulation_time, return mse end +# Calculate the mean relative error (MRE) between interpolated simulation values and reference values +# over a common time range. +# +# Arguments: +# - `reference_time` : Time points for the reference data. +# - `reference_values` : Data points for the reference data. +# - `simulation_time` : Time points for the simulation data. +# - `simulation_values` : Data points for the simulation data. +function interpolated_mre(reference_time, reference_values, simulation_time, simulation_values) + if last(simulation_time) > last(reference_time) + @warn "simulation time range is larger than reference time range. " * + "Only checking values within reference time range." + end + + # Remove reference time points outside the simulation time + start = searchsortedfirst(reference_time, first(simulation_time)) + end_ = searchsortedlast(reference_time, last(simulation_time)) + common_time_range = reference_time[start:end_] + + # Interpolate simulation data at the common time points + interpolated_values = [linear_interpolation(simulation_time, simulation_values, t) + for t in common_time_range] + + filtered_values = reference_values[start:end_] + + # Calculate MRE only over the common time range + mre = sum(abs.(interpolated_values .- filtered_values) ./ abs.(filtered_values)) / length(common_time_range) + return mre +end + + function extract_number_from_filename(filename) # This regex matches the last sequence of digits in the filename m = match(r"(\d+)(?!.*\d)", filename)