Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@svchb
svchb committed Jul 17, 2024
1 parent 61e3325 commit c53246d
Show file tree
Hide file tree
Showing 8 changed files with 136 additions and 29 deletions.
11 changes: 9 additions & 2 deletions examples/fluid/falling_water_spheres_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,13 @@ sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calc
surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
correction=AkinciFreeSurfaceCorrection(fluid_density))

# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator,
# state_equation, fluid_smoothing_kernel,
# fluid_smoothing_length,
# viscosity=viscosity,
# acceleration=(0.0, -gravity),
# correction=AkinciFreeSurfaceCorrection(fluid_density))

# sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
# fluid_smoothing_length,
# sound_speed, viscosity=viscosity,
Expand Down Expand Up @@ -85,11 +92,11 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
viscosity=ViscosityAdami(nu=wall_viscosity))

boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
adhesion_coefficient=0.001)
adhesion_coefficient=0.1)

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

info_callback = InfoCallback(interval=50)
Expand Down
12 changes: 5 additions & 7 deletions examples/fluid/test copy.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,6 @@ saving_callback = SolutionSavingCallback(dt=0.005, output_directory="out",

stepsize_callback = StepsizeCallback(cfl=2.5)


callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback)

# Use a Runge-Kutta method with automatic (error based) time step size control.
Expand All @@ -193,16 +192,15 @@ sol = solve(ode, CKLLSRK54_3M_3R(),
# sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
# dt=1.0, # This is overwritten by the stepsize callback
# save_everystep=false, callback=callbacks);
#~157-179s to 1%
# cfl 2.0 diverged
# cfl 1.5 19700s to 100% 750-780s to 5%

#~157-179s to 1%
# cfl 2.0 diverged
# cfl 1.5 19700s to 100% 750-780s to 5%

# sol = solve(ode, Tsit5(),
# dt=1.0, # This is overwritten by the stepsize callback
# save_everystep=false, callback=callbacks);
#~219-246s to 1%
# cfl 2.0 842-850s to 5%
#~219-246s to 1%
# cfl 2.0 842-850s to 5%

# sol = solve(ode, ParsaniKetchesonDeconinck3S32(),
# dt=1.0, # This is overwritten by the stepsize callback
Expand Down
5 changes: 2 additions & 3 deletions examples/fluid/test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,15 +75,14 @@ viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
# 0.011),
# correction=AkinciFreeSurfaceCorrection(fluid_density))


sphere_surface_tension = EntropicallyDampedSPHSystem(water, fluid_smoothing_kernel,
fluid_smoothing_length,
sound_speed, viscosity=viscosity,
density_calculator=ContinuityDensity(),
acceleration=(0.0, -gravity),
surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.005),
surface_normal_method=AkinciSurfaceNormal(smoothing_kernel=WendlandC6Kernel{2}(),
smoothing_length=4 *
surface_normal_method=AkinciSurfaceNormal(smoothing_kernel=WendlandC6Kernel{2}(),
smoothing_length=4 *
fluid_particle_spacing))

# sphere_surface_tension = WeaklyCompressibleSPHSystem(water, fluid_density_calculator,
Expand Down
3 changes: 3 additions & 0 deletions src/general/semidiscretization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,6 +277,9 @@ function semidiscretize(semi, tspan; reset_threads=true, data_type=nothing)
# Get the neighborhood search for this system
neighborhood_search = get_neighborhood_search(system, semi)

PointNeighbors.initialize!(neighborhood_search, initial_coordinates(system),
initial_coordinates(system))

# Initialize this system
initialize!(system, neighborhood_search)

Expand Down
3 changes: 2 additions & 1 deletion src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,

cache = (; create_cache_model(viscosity, n_particles, NDIMS)...,
create_cache_model(initial_density, density_calculator)...,
create_cache_model(correction, initial_density, NDIMS, n_particles)...)
create_cache_model(correction, initial_density, NDIMS, n_particles)...,
(; colorfield=zeros(eltype(smoothing_length), n_particles))...)

return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation,
density_calculator, smoothing_kernel,
Expand Down
14 changes: 14 additions & 0 deletions src/schemes/boundary/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -390,3 +390,17 @@ end
function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem)
return Inf
end

function initialize!(system::BoundarySPHSystem, neighborhood_search)
# TODO: dispatch on boundary model
system_coords = system.coordinates
(; smoothing_kernel, smoothing_length) = system.boundary_model

foreach_point_neighbor(system, system,
system_coords, system_coords,
neighborhood_search, points=eachparticle(system)) do particle, neighbor, pos_diff, distance
system.boundary_model.cache.colorfield[particle] += kernel(smoothing_kernel,
distance,
smoothing_length)
end
end
114 changes: 99 additions & 15 deletions src/schemes/fluid/surface_normal_sph.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,29 +10,27 @@ 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,
surface_tension::SurfaceTensionAkinci, u_system,
v_neighbor_system, u_neighbor_system,
semi, surfn)
function calc_normal_akinci!(system, 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)

# println("sml ", smoothing_length != system.smoothing_length)
# println("smk ", smoothing_kernel !== system.smoothing_kernel)
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))
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
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)
Expand All @@ -50,9 +48,95 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem,
return system
end

function calc_normal_akinci!(system, neighbor_system, surface_tension, u_system,
v_neighbor_system, u_neighbor_system,
semi, surfn)
# 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)
(; 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)
# nhs_bnd = get_neighborhood_search(neighbor_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))
# nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius,
# nparticles(neighbor_system))
# PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords)
# PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords)
# end

# First we need to calculate the smoothed colorfield values
# TODO: move colorfield to extra step
# TODO: this is only correct for a single fluid
colorfield = zeros(eltype(neighbor_system), nparticles(neighbor_system))
# colorfield_bnd = zeros(eltype(neighbor_system), nparticles(neighbor_system))

foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
colorfield[neighbor] += kernel(smoothing_kernel, distance, smoothing_length)
end

# foreach_point_neighbor(neighbor_system, neighbor_system,
# neighbor_system_coords, neighbor_system_coords,
# nhs_bnd) do particle, neighbor, pos_diff, distance
# println("test")
# colorfield_bnd[particle] += kernel(smoothing_kernel, distance, smoothing_length)
# end

# Since we don't want to calculate the unused boundary weight sum we normalize against the maximum value
colorfield = colorfield ./ (2.0 * maximum(colorfield))

# foreach_point_neighbor(system, neighbor_system,
# system_coords, neighbor_system_coords,
# nhs) do particle, neighbor, pos_diff, distance
# colorfield[neighbor] = colorfield[neighbor] /
# (colorfield[neighbor] +
# neighbor_system.boundary_model.cache.colorfield[neighbor])
# end

# colorfield = colorfield / (colorfield + neighbor_system.boundary_model.cache.colorfield)

# for i in 1:nparticles(neighbor_system)
# if colorfield[i] > eps()
# println(colorfield)
# break
# end
# end

foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
if colorfield[neighbor] > 0.05
m_b = hydrodynamic_mass(system, particle)
density_neighbor = particle_density(v, system, particle)
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
cache.surface_normal[i, particle] = 0.0
# println(m_b / density_neighbor * grad_kernel[i] * smoothing_length)
end
cache.neighbor_count[particle] += 1
end
end

return system
end

function calc_normal_akinci!(system, neighbor_system, u_system, v, v_neighbor_system,
u_neighbor_system, semi, surfn)
# Normal not needed
return system
end
Expand Down Expand Up @@ -80,8 +164,8 @@ function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, se
return system
end

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

# Reset surface normal
Expand All @@ -92,7 +176,7 @@ function compute_surface_normal!(system, surface_tension::SurfaceTensionAkinci,
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, surface_tension, u, v_neighbor_system,
calc_normal_akinci!(system, neighbor_system, u, v, v_neighbor_system,
u_neighbor_system, semi, surface_normal_method)
remove_invalid_normals!(system, surface_tension)
end
Expand Down
3 changes: 2 additions & 1 deletion src/schemes/fluid/weakly_compressible_sph/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,8 @@ function WeaklyCompressibleSPHSystem(initial_condition,
acceleration_, viscosity,
density_diffusion, correction,
pressure_acceleration,
source_terms, surface_tension, surface_normal_method, buffer, cache)
source_terms, surface_tension, surface_normal_method,
buffer, cache)
end

create_cache_wcsph(correction, density, NDIMS, nparticles) = (;)
Expand Down

0 comments on commit c53246d

Please sign in to comment.