From c53246dbca9f08820e61ad736794c591697b6baf Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 17 Jul 2024 17:33:59 +0200 Subject: [PATCH] update --- examples/fluid/falling_water_spheres_2d.jl | 11 +- examples/fluid/test copy.jl | 12 +- examples/fluid/test.jl | 5 +- src/general/semidiscretization.jl | 3 + .../dummy_particles/dummy_particles.jl | 3 +- src/schemes/boundary/system.jl | 14 +++ src/schemes/fluid/surface_normal_sph.jl | 114 +++++++++++++++--- .../fluid/weakly_compressible_sph/system.jl | 3 +- 8 files changed, 136 insertions(+), 29 deletions(-) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 1d1565809..0ac0f18ca 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -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, @@ -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) diff --git a/examples/fluid/test copy.jl b/examples/fluid/test copy.jl index 6ada96f64..b0a3a28aa 100644 --- a/examples/fluid/test copy.jl +++ b/examples/fluid/test copy.jl @@ -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. @@ -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 diff --git a/examples/fluid/test.jl b/examples/fluid/test.jl index c44fe990b..0b60817f1 100644 --- a/examples/fluid/test.jl +++ b/examples/fluid/test.jl @@ -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, diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index b1313dd44..60ea6a4b7 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -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) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 6b544aa78..04e7f0a99 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -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, diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 8f181ce7e..208e88909 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -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 diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 2ee1eb75f..810fd0180 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -10,10 +10,8 @@ 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 @@ -21,18 +19,18 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, 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) @@ -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 @@ -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 @@ -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 diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index e82ab54b8..340cc0173 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -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) = (;)