Skip to content

Commit

Permalink
inter
Browse files Browse the repository at this point in the history
  • Loading branch information
@svchb
svchb committed Aug 1, 2024
1 parent 1969525 commit cf28d2e
Show file tree
Hide file tree
Showing 9 changed files with 275 additions and 54 deletions.
94 changes: 94 additions & 0 deletions examples/fluid/falling_water_spheres_morris_2d.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
# In this example two circles of water drop to the floor demonstrating the difference
# between the behavior with and without surface tension modelling.
using TrixiParticles
using OrdinaryDiffEq

# ==========================================================================================
# ==== Resolution
fluid_particle_spacing = 0.005

boundary_layers = 3
spacing_ratio = 1

# ==========================================================================================
# ==== Experiment Setup
gravity = 9.81
tspan = (0.0, 0.3)

# Boundary geometry and initial fluid particle positions
initial_fluid_size = (0.0, 0.0)
tank_size = (2.0, 0.5)

fluid_density = 1000.0
sound_speed = 100
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)

sphere_radius = 0.05

sphere1_center = (0.5, 0.2)
sphere2_center = (1.5, 0.2)
sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center,
fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0))
sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,
fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0))

# ==========================================================================================
# ==== Fluid
fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()

fluid_density_calculator = ContinuityDensity()

nu = 0.005
alpha = 8 * nu / (fluid_smoothing_length * sound_speed)
viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1)

sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
fluid_smoothing_length,
sound_speed, viscosity=viscosity,
density_calculator=ContinuityDensity(),
acceleration=(0.0, -gravity),
surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.05))

sphere = WeaklyCompressibleSPHSystem(sphere2, 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()
wall_viscosity = nu
boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
state_equation=state_equation,
boundary_density_calculator,
fluid_smoothing_kernel, fluid_smoothing_length,
viscosity=ViscosityAdami(nu=wall_viscosity))

boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
adhesion_coefficient=1.0)

# ==========================================================================================
# ==== Simulation
semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system)
ode = semidiscretize(semi, tspan)

info_callback = InfoCallback(interval=50)
saving_callback = SolutionSavingCallback(dt=0.01, 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, RDPK3SpFSAL35(),
abstol=1e-7, # Default abstol is 1e-6
reltol=1e-4, # Default reltol is 1e-3
save_everystep=false, callback=callbacks);
2 changes: 1 addition & 1 deletion src/TrixiParticles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure,
max_density, min_density, avg_density
export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d,
interpolate_plane_2d_vtk
export SurfaceTensionAkinci, CohesionForceAkinci
export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris
export ColorfieldSurfaceNormal

end # module
3 changes: 2 additions & 1 deletion src/schemes/fluid/entropically_damped_sph/rhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ function interact!(dv, v_particle_system, u_particle_system,

dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b,
particle_system, neighbor_system,
particle, neighbor, pos_diff, distance)
particle, neighbor, pos_diff, distance,
rho_a, rho_b)

dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
particle, neighbor, pos_diff, distance)
Expand Down
19 changes: 7 additions & 12 deletions src/schemes/fluid/entropically_damped_sph/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,10 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V,

cache = create_cache_density(initial_condition, density_calculator)
cache = (;
create_cache_edac(surface_tension, ELTYPE, NDIMS, n_particles)...,
create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS,
n_particles)...,
create_cache_surface_tension(surface_tension, ELTYPE, NDIMS,
n_particles)...,
cache...)

new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass),
Expand All @@ -119,16 +122,6 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V,
end
end

function create_cache_edac(surface_tension_model, ELTYPE, NDIMS, nparticles)
return (;)
end

function create_cache_edac(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles)
surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles)
neighbor_count = Array{ELTYPE, 1}(undef, nparticles)
return (; surface_normal, neighbor_count)
end

function Base.show(io::IO, system::EntropicallyDampedSPHSystem)
@nospecialize system # reduce precompilation time

Expand Down Expand Up @@ -205,7 +198,9 @@ end
function update_quantities!(system::EntropicallyDampedSPHSystem, v, u,
v_ode, u_ode, semi, t)
compute_density!(system, u, u_ode, semi, system.density_calculator)
compute_surface_normal!(system, system.surface_tension, v, u, v_ode, u_ode, semi, t)
compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi,
t)
compute_curvature!(system, system.surface_tension, v, u, v_ode, u_ode, semi, t)
end

function write_v0!(v0, system::EntropicallyDampedSPHSystem, density_calculator)
Expand Down
131 changes: 110 additions & 21 deletions src/schemes/fluid/surface_normal_sph.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,26 @@ function ColorfieldSurfaceNormal(; smoothing_kernel, smoothing_length)
return ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length)
end

function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles)
return (;)
end

function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles)
surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles)
neighbor_count = Array{ELTYPE, 1}(undef, nparticles)
return (; surface_normal, neighbor_count)
end

function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system,
u_neighbor_system, semi, surfn)
# Normal not needed
return system
end

# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids"
# Note: This is the simplest form of normal approximation commonly used in SPH and comes
# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries
function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v,
v_neighbor_system, u_neighbor_system, semi, surfn)
# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics"
function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v,
v_neighbor_system, u_neighbor_system, semi, surfn)
(; cache) = system
(; smoothing_kernel, smoothing_length) = surfn

Expand All @@ -37,9 +52,8 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v,
grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length)
for i in 1:ndims(system)
# The `smoothing_length` here is used for scaling
# TODO move this to the surface tension model since this is not a general thing
cache.surface_normal[i, particle] += m_b / density_neighbor *
grad_kernel[i] * smoothing_length
grad_kernel[i]
end

cache.neighbor_count[particle] += 1
Expand All @@ -49,10 +63,9 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v,
end

# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids"
# Note: This is the simplest form of normal approximation commonly used in SPH and comes
# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries
function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system, v,
v_neighbor_system, u_neighbor_system, semi, surfn)
# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics"
function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_system, v,
v_neighbor_system, u_neighbor_system, semi, surfn)
(; cache) = system
(; colorfield, colorfield_bnd) = neighbor_system.boundary_model.cache
(; smoothing_kernel, smoothing_length) = surfn
Expand Down Expand Up @@ -101,9 +114,8 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system,
smoothing_length)
for i in 1:ndims(system)
# The `smoothing_length` here is used for scaling
# TODO move this to the surface tension model since this is not a general thing
cache.surface_normal[i, particle] += m_b / density_neighbor *
grad_kernel[i] * smoothing_length
grad_kernel[i]
end
cache.neighbor_count[particle] += 1
end
Expand All @@ -112,8 +124,7 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system,
return system
end

function calc_normal_akinci!(system, neighbor_system, u_system, v, v_neighbor_system,
u_neighbor_system, semi, surfn)
function remove_invalid_normals!(system, surface_tension)
# Normal not needed
return system
end
Expand All @@ -132,18 +143,36 @@ function remove_invalid_normals!(system::FluidSystem, surface_tension::SurfaceTe
return system
end

function remove_invalid_normals!(system, surface_tension)
# Normal not needed
function remove_invalid_normals!(system::FluidSystem, surface_tension::SurfaceTensionMorris)
(; cache, smoothing_length) = system

# We remove invalid normals i.e. they have a small norm (eq. 20)
normal_condition2 = (0.01 / smoothing_length)^2

for particle in each_moving_particle(system)
particle_surface_normal = cache.surface_normal[1:ndims(system), particle]
norm2 = dot(particle_surface_normal, particle_surface_normal)

# see eq. 21
if norm2 > normal_condition2
cache.surface_normal[1:ndims(system), particle] = particle_surface_normal /
sqrt(norm2)
else
cache.surface_normal[1:ndims(system), particle] .= 0
end
end

return system
end

function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t)
return system
end

function compute_surface_normal!(system::FluidSystem, surface_tension::SurfaceTensionAkinci,
function compute_surface_normal!(system::FluidSystem,
surface_normal_method::ColorfieldSurfaceNormal,
v, u, v_ode, u_ode, semi, t)
(; cache, surface_normal_method) = system
(; cache, surface_tension) = system

# Reset surface normal
set_zero!(cache.surface_normal)
Expand All @@ -153,9 +182,69 @@ function compute_surface_normal!(system::FluidSystem, surface_tension::SurfaceTe
u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)

calc_normal_akinci!(system, neighbor_system, u, v, v_neighbor_system,
u_neighbor_system, semi, surface_normal_method)
calc_normal!(system, neighbor_system, u, v, v_neighbor_system,
u_neighbor_system, semi, surface_normal_method)
remove_invalid_normals!(system, surface_tension)
end
return system
end

# Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics"
function calc_curvature!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v,
v_neighbor_system, u_neighbor_system, semi, surfn)
(; cache) = system
(; smoothing_kernel, smoothing_length) = surfn

system_coords = current_coordinates(u_system, system)
neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system)
nhs = get_neighborhood_search(system, neighbor_system, semi)

if smoothing_length != system.smoothing_length ||
smoothing_kernel !== system.smoothing_kernel
# TODO: this is really slow but there is no way to easily implement multiple search radia
search_radius = compact_support(smoothing_kernel, smoothing_length)
nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius,
nparticles(system))
PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords)
end

foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
m_b = hydrodynamic_mass(neighbor_system, neighbor)
density_neighbor = particle_density(v_neighbor_system,
neighbor_system, neighbor)
grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length)
for i in 1:ndims(system)
# The `smoothing_length` here is used for scaling
cache.surface_normal[i, particle] += m_b / density_neighbor *
grad_kernel[i]
end

cache.neighbor_count[particle] += 1
end

return system
end

function compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
return system
end

function compute_curvature!(system::FluidSystem,
surface_tension::ViscosityMorris,
v, u, v_ode, u_ode, semi, t)
(; cache, surface_tension) = system

# Reset surface curvature
set_zero!(cache.curvature)

@trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system
u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)

calc_curvature!(system, neighbor_system, u, v, v_neighbor_system,
u_neighbor_system, semi, surface_normal_method)
end
return system
end
Loading

0 comments on commit cf28d2e

Please sign in to comment.