diff --git a/examples/fluid/particle_refinement_2d.jl b/examples/fluid/particle_refinement_2d.jl new file mode 100644 index 000000000..d873d43f0 --- /dev/null +++ b/examples/fluid/particle_refinement_2d.jl @@ -0,0 +1,86 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (2.0, 2.0) +tank_size = (2.0, 2.0) + +fluid_density = 1000.0 +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=2.0) + +sphere = SphereShape(0.5 * fluid_particle_spacing, 0.3, (1.0, 1.0), + fluid_density, sphere_type=RoundSphere()) +# ========================================================================================== +# ==== Fluid + +fluid = setdiff(tank.fluid, sphere) + +smoothing_length = 1.2 * fluid_particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() +pressure = 1000.0 +particle_refinement = ParticleRefinement(; + refinement_pattern=TrixiParticles.HexagonalSplitting(), + max_spacing_ratio=1.2) +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed; viscosity=ViscosityAdami(; nu=1e-4), + particle_refinement=particle_refinement, + transport_velocity=TransportVelocityAdami(pressure), + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary + +# This is to set another boundary density calculation with `trixi_include` +boundary_density_calculator = AdamiPressureExtrapolation() + +# This is to set wall viscosity with `trixi_include` +viscosity_wall = nothing +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing) + +boundary_model_sphere = BoundaryModelDummyParticles(sphere.density, sphere.mass, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) +boundary_system_sphere = BoundarySPHSystem(sphere, boundary_model_sphere, movement=nothing) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system, boundary_system_sphere) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +# This is to easily add a new callback with `trixi_include` +extra_callback = nothing + +refinement_callback = ParticleRefinementCallback(interval=1) + +callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback(), + refinement_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c316cd801..b84d9fdfb 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -22,7 +22,8 @@ using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, - get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem + get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, + RecursiveArrayTools @reexport using StaticArrays: SVector using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt @@ -45,6 +46,7 @@ include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") include("schemes/schemes.jl") +include("refinement/refinement.jl") # Note that `semidiscretization.jl` depends on the system types and has to be # included separately. `gpu.jl` in turn depends on the semidiscretization type. @@ -52,6 +54,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("refinement/resize.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition @@ -59,7 +62,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, OutFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, ParticleRefinementCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, @@ -72,7 +75,7 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BernoulliPressureExtrapolation - +export ParticleRefinement, SpatialRefinementCriterion export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 1aefdff72..3d8f5701a 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") include("update.jl") +include("refinement.jl") diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl new file mode 100644 index 000000000..8f04bed40 --- /dev/null +++ b/src/callbacks/refinement.jl @@ -0,0 +1,126 @@ +struct ParticleRefinementCallback{I} + interval::I +end + +function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) + if dt > 0 && interval !== -1 + throw(ArgumentError("setting both `interval` and `dt` is not supported")) + end + + # Update in intervals in terms of simulation time + if dt > 0 + interval = Float64(dt) + + # Update every time step (default) + elseif interval == -1 + interval = 1 + end + + refinement_callback = ParticleRefinementCallback(interval) + + if dt > 0 + # Add a `tstop` every `dt`, and save the final solution. + return PeriodicCallback(refinement_callback, dt, + initialize=initial_refinement!, + save_positions=(false, false)) + else + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(refinement_callback, refinement_callback, + initialize=initial_refinement!, + save_positions=(false, false)) + end +end + +# initialize +function initial_refinement!(cb, u, t, integrator) + # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`) + # or `cb.affect!.affect!` (with `PeriodicCallback`). + # Let recursive dispatch handle this. + + initial_refinement!(cb.affect!, u, t, integrator) +end + +function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator) + semi = integrator.p + + foreach_system(semi) do system + resize_refinement!(system) + end + + cb(integrator) +end + +# condition +function (refinement_callback::ParticleRefinementCallback)(u, t, integrator) + (; interval) = refinement_callback + + return condition_integrator_interval(integrator, interval) +end + +# affect +function (refinement_callback::ParticleRefinementCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # Basically `get_tmp_cache(integrator)` to write into in order to be non-allocating + # https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Caches + v_tmp, u_tmp = integrator.cache.tmp.x + + v_tmp .= v_ode + u_tmp .= u_ode + + refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + + resize!(integrator, (length(v_ode), length(u_ode))) + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +Base.resize!(a::RecursiveArrayTools.ArrayPartition, sizes::Tuple) = resize!.(a.x, sizes) + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")") +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(dt=", cb.affect!.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect! + setup = [ + "interval" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect!.affect! + setup = [ + "dt" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end diff --git a/src/general/corrections.jl b/src/general/corrections.jl index b1c66a06a..c63c25a27 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -220,7 +220,7 @@ function compute_correction_values!(system, kernel_correction_coefficient[particle] += volume * smoothing_kernel(system, distance) if distance > sqrt(eps()) - tmp = volume * smoothing_kernel_grad(system, pos_diff, distance) + tmp = volume * smoothing_kernel_grad(system, pos_diff, distance, particle) for i in axes(dw_gamma, 1) dw_gamma[i, particle] += tmp[i] end @@ -311,7 +311,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, pos_diff, distance volume = mass[neighbor] / density_fun(neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) iszero(grad_kernel) && return @@ -349,7 +349,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, smoothing_length, system, particle) - return smoothing_kernel_grad(system, pos_diff, distance) + return smoothing_kernel_grad(system, pos_diff, distance, particle) end # Compute gradient of corrected kernel diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 2a3d45ea0..e15eb33e0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -66,7 +66,7 @@ function summation_density!(system, semi, u, u_ode, density; points=particles) do particle, neighbor, pos_diff, distance mass = hydrodynamic_mass(neighbor_system, neighbor) - density[particle] += mass * smoothing_kernel(system, distance) + density[particle] += mass * smoothing_kernel(system, distance, particle) end end end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..f734a6191 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -73,14 +73,7 @@ function Semidiscretization(systems...; # Other checks might be added here later. check_configuration(systems) - sizes_u = [u_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) - for i in eachindex(sizes_u)) - sizes_v = [v_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) - for i in eachindex(sizes_v)) + ranges_v, ranges_u = ranges_vu(systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. @@ -92,6 +85,16 @@ function Semidiscretization(systems...; return Semidiscretization(systems, ranges_u, ranges_v, searches) end +function ranges_vu(systems) + sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_u = [(sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) for i in eachindex(sizes_u)] + + sizes_v = [v_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_v = [(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)] + + return ranges_v, ranges_u +end + # Inline show function e.g. Semidiscretization(neighborhood_search=...) function Base.show(io::IO, semi::Semidiscretization) @nospecialize semi # reduce precompilation time @@ -134,8 +137,8 @@ function create_neighborhood_search(neighborhood_search, system, neighbor) end @inline function compact_support(system, neighbor) - (; smoothing_kernel, smoothing_length) = system - return compact_support(smoothing_kernel, smoothing_length) + (; smoothing_kernel) = system + return compact_support(smoothing_kernel, smoothing_length(system, 1)) # TODO end @inline function compact_support(system::OpenBoundarySPHSystem, neighbor) @@ -407,7 +410,9 @@ end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) (; systems) = semi - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) + return 1.0 + # TODO + # return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) end function drift!(du_ode, v_ode, u_ode, semi, t) diff --git a/src/general/system.jl b/src/general/system.jl index 48d556eba..88b9e7bf2 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -106,30 +106,23 @@ end @inline set_particle_pressure!(v, system, particle, pressure) = v -@inline function smoothing_kernel(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel(smoothing_kernel, distance, smoothing_length) +@inline function smoothing_kernel(system, distance, particle) + (; smoothing_kernel) = system + return kernel(smoothing_kernel, distance, smoothing_length(system, particle)) end -@inline function smoothing_kernel_deriv(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel_deriv(smoothing_kernel, distance, smoothing_length) -end - -@inline function smoothing_kernel_grad(system, pos_diff, distance) - return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) -end - -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance) +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) (; smoothing_kernel, smoothing_length) = system.boundary_model return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) end @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) - return corrected_kernel_grad(system.smoothing_kernel, pos_diff, distance, - system.smoothing_length, system.correction, system, - particle) + (; smoothing_kernel) = system + + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length(system, particle), + system.correction, system, particle) end # System update orders. This can be dispatched if needed. diff --git a/src/refinement/merge.jl b/src/refinement/merge.jl new file mode 100644 index 000000000..ea55bb0bc --- /dev/null +++ b/src/refinement/merge.jl @@ -0,0 +1,121 @@ +function merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + merge_particles!(system, semi, v, u) + end + + deleteat!(semi, v_ode, u_ode, v_tmp, u_tmp) + + return semi +end + +@inline merge_particles!(system, semi, v, u) = system + +@inline function merge_particles!(system::FluidSystem, semi, v, u) + return merge_particles!(system, system.particle_refinement, semi, v, u) +end + +@inline merge_particles!(system::FluidSystem, ::Nothing, semi, v, u) = system + +@inline function merge_particles!(system::FluidSystem, particle_refinement, semi, v, u) + (; delete_candidates) = system.particle_refinement + + delete_candidates .= false + + # Merge particles iteratively + for _ in 1:3 + merge_particles_inner!(system, particle_refinement, semi, v, u) + end + + return system +end + +function merge_particles_inner!(system, particle_refinement, semi, v, u) + (; smoothing_kernel, cache) = system + (; mass_ref, max_spacing_ratio, merge_candidates, delete_candidates) = particle_refinement + + set_zero!(merge_candidates) + + system_coords = current_coordinates(u, system) + neighborhood_search = get_neighborhood_search(system, semi) + + # Collect merge candidates + foreach_point_neighbor(system, system, system_coords, system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + particle == neighbor && return + delete_candidates[particle] && return + + m_a = hydrodynamic_mass(system, particle) + m_b = hydrodynamic_mass(system, neighbor) + + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a <= m_max + m_merge = m_a + m_b + m_max_min = min(m_max, max_spacing_ratio * mass_ref[neighbor]) + if m_merge < m_max_min + if merge_candidates[particle] == 0 + merge_candidates[particle] = neighbor + else + stored_neighbor = current_coords(u, system, merge_candidates[particle]) + pos_diff_stored = stored_neighbor - current_coords(u, system, particle) + + if distance < norm(pos_diff_stored) + merge_candidates[particle] = neighbor + end + end + end + end + end + + # Merge and delete particles + for particle in findall(!iszero, merge_candidates) + + candidate = merge_candidates[particle] + + if particle == merge_candidates[candidate] + if particle < candidate + m_a = hydrodynamic_mass(system, particle) + m_b = hydrodynamic_mass(system, candidate) + + m_merge = m_a + m_b + + pos_a = current_coords(u, system, particle) + pos_b = current_coords(u, system, candidate) + + vel_a = current_velocity(v, system, particle) + vel_b = current_velocity(v, system, candidate) + + pos_merge = (m_a * pos_a + m_b * pos_b) / m_merge + vel_merge = (m_a * vel_a + m_b * vel_b) / m_merge + + for dim in 1:ndims(system) + u[dim, particle] = pos_merge[dim] + v[dim, particle] = vel_merge[dim] + end + + # update smoothing length + h_a = smoothing_length(system, particle) + h_b = smoothing_length(system, candidate) + + # TODO: + # Check normalization_factor = smoothing_kernel(zero(pos_diff), one(h_a)) + tmp_m = m_merge * normalization_factor(smoothing_kernel, h_a) + + tmp_a = m_a * kernel(smoothing_kernel, norm(pos_merge - pos_a), h_a) + tmp_b = m_b * kernel(smoothing_kernel, norm(pos_merge - pos_b), h_b) + + cache.smoothing_length[particle] = (tmp_m / (tmp_a + tmp_b))^(1 / + ndims(system)) + + system.mass[particle] = m_merge + else + delete_candidates[candidate] = true + end + end + end + + return system +end diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl new file mode 100644 index 000000000..399a42bab --- /dev/null +++ b/src/refinement/refinement.jl @@ -0,0 +1,207 @@ +include("refinement_criteria.jl") +include("refinement_pattern.jl") +include("split.jl") +include("merge.jl") + +struct ParticleRefinement{SP, RC, ELTYPE} + refinement_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles + split_candidates :: Vector{Int} # variable length + n_new_particles :: Ref{Int} +end + +function ParticleRefinement(; refinement_pattern, max_spacing_ratio, + refinement_criteria=SpatialRefinementCriterion()) + mass_ref = Vector{eltype(max_spacing_ratio)}() + delete_candidates = Vector{Bool}() + + if !(refinement_criteria isa Tuple) + refinement_criteria = (refinement_criteria,) + end + + return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, + mass_ref, Int[], delete_candidates, Int[], Ref(0)) +end + +resize_refinement!(system) = system + +function resize_refinement!(system::FluidSystem) + resize_refinement!(system, system.particle_refinement) +end + +resize_refinement!(system, ::Nothing) = system + +function resize_refinement!(system, particle_refinement) + resize!(particle_refinement.mass_ref, nparticles(system)) + resize!(particle_refinement.delete_candidates, nparticles(system)) + resize!(particle_refinement.merge_candidates, nparticles(system)) + + particle_refinement.mass_ref .= system.mass + + return system +end + +function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + # Check refinement criteria + check_refinement_criteria!(semi, v_ode, u_ode) + + # Update spacing of particles (Algorthm 1) + update_particle_spacing(semi, v_ode, u_ode) + + # Split particles (Algorithm 2) + split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) + + # Merge particles (Algorithm 3) + merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) + + # Shift particles + # TODO + + # Correct particles + # TODO + + # Update smoothing lengths + upate_smoothing_lengths!(semi, u_ode) + + # Resize neighborhood search + # TODO + + return semi +end + +function update_particle_spacing(semi, v_ode, u_ode) + foreach_system(semi) do system + update_particle_spacing(system, v_ode, u_ode, semi) + end +end + +# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl` +@inline update_particle_spacing(system, v_ode, u_ode, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, v_ode, u_ode, semi) + update_particle_spacing(system, system.particle_refinement, v_ode, u_ode, semi) +end + +@inline update_particle_spacing(system::FluidSystem, ::Nothing, v_ode, u_ode, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, + v_ode, u_ode, semi) + (; smoothing_length, smoothing_length_factor) = system.cache + (; mass_ref, max_spacing_ratio) = particle_refinement + + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + + system_coords = current_coordinates(u, system) + + for particle in eachparticle(system) + dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, u_ode, system_coords, + particle) + + if dp_max / dp_min < max_spacing_ratio^3 + new_spacing = min(dp_max, max_spacing_ratio * dp_min) + else + new_spacing = dp_avg + end + + smoothing_length[particle] = smoothing_length_factor * new_spacing + mass_ref[particle] = particle_density(v, system, particle) * + new_spacing^(ndims(system)) + end + + return system +end + +function upate_smoothing_lengths!(semi, u_ode) + foreach_system(semi) do system + u = wrap_u(u_ode, system, semi) + upate_smoothing_lengths!(system, semi, u) + end + + return semi +end + +upate_smoothing_lengths!(system, semi, u) = system + +function upate_smoothing_lengths!(system::FluidSystem, semi, u) + upate_smoothing_lengths!(system, system.particle_refinement, semi, u) +end + +upate_smoothing_lengths!(system::FluidSystem, ::Nothing, semi, u) = system + +function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) + (; smoothing_length, smoothing_length_factor) = system.cache + + neighborhood_search = get_neighborhood_search(system, semi) + + system_coords = current_coordinates(u, system) + + # TODO: Think about a solution when density is not constant in `initial_condition` + # Problem: Resized systems + rho_0 = first(system.initial_condition.density) + + for particle in each_moving_particle(system) + mass_avg = zero(eltype(system)) + counter_neighbors = 0 + + PointNeighbors.foreach_neighbor(system_coords, system_coords, neighborhood_search, + particle) do particle, neighbor, pos_diff, distance + mass_avg += hydrodynamic_mass(system, neighbor) + + counter_neighbors += 1 + end + + mass_avg /= counter_neighbors + + smoothing_length[particle] = smoothing_length_factor * + (mass_avg / rho_0)^(1 / ndims(system)) + end + + return system +end + +@inline function min_max_avg_spacing(system, semi, u_ode, system_coords, particle) + dp_min = particle_spacing(system, particle) + dp_max = particle_spacing(system, particle) + dp_avg = particle_spacing(system, particle) + counter_neighbors = 0 + + foreach_system(semi) do neighbor_system + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) + + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + PointNeighbors.foreach_neighbor(system_coords, neighbor_coords, neighborhood_search, + particle) do particle, neighbor, pos_diff, distance + dp_neighbor = particle_spacing(neighbor_system, neighbor) + + dp_min = min(dp_min, dp_neighbor) + dp_max = max(dp_max, dp_neighbor) + dp_avg += dp_neighbor + + counter_neighbors += 1 + end + end + + dp_avg = counter_neighbors == 0 ? dp_avg : dp_avg / counter_neighbors + + return dp_min, dp_max, dp_avg +end + +@inline particle_spacing(system, particle) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system::FluidSystem, particle) + return particle_spacing(system, system.particle_refinement, particle) +end + +@inline particle_spacing(system, ::Nothing, _) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system, refinement, particle) + (; smoothing_length_factor) = system.cache + return smoothing_length(system, particle) / smoothing_length_factor +end diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl new file mode 100644 index 000000000..b8106cbb8 --- /dev/null +++ b/src/refinement/refinement_criteria.jl @@ -0,0 +1,69 @@ +abstract type RefinementCriteria end + +struct SpatialRefinementCriterion <: RefinementCriteria end + +struct SolutionRefinementCriterion <: RefinementCriteria end + +function check_refinement_criteria!(semi, v_ode, u_ode) + foreach_system(semi) do system + check_refinement_criteria!(system, v_ode, u_ode, semi) + end +end + +@inline check_refinement_criteria!(system, v_ode, u_ode, semi) = system + +@inline function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + check_refinement_criteria!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) +end + +@inline function check_refinement_criteria!(system::FluidSystem, ::Nothing, + v, u, v_ode, u_ode, semi) + system +end + +@inline function check_refinement_criteria!(system::FluidSystem, refinement, + v, u, v_ode, u_ode, semi) + (; refinement_criteria) = system.particle_refinement + for criterion in refinement_criteria + criterion(system, v, u, v_ode, u_ode, semi) + end +end + +@inline function (criterion::SpatialRefinementCriterion)(system, v, u, v_ode, u_ode, semi) + system_coords = current_coordinates(u, system) + + foreach_system(semi) do neighbor_system + set_particle_spacing!(system, neighbor_system, system_coords, v_ode, u_ode, semi) + end + return system +end + +@inline set_particle_spacing!(system, _, _, _, _, _) = system + +@inline function set_particle_spacing!(particle_system, + neighbor_system::Union{BoundarySystem, SolidSystem}, + system_coords, v_ode, u_ode, semi) + (; smoothing_length, smoothing_length_factor) = particle_system.cache + + neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, semi) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_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 + + dp_particle = particle_spacing(particle_system, particle) + dp_neighbor = particle_spacing(neighbor_system, neighbor) + + smoothing_length[particle] = smoothing_length_factor * min(dp_neighbor, dp_particle) + end + + return particle_system +end diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl new file mode 100644 index 000000000..c01022314 --- /dev/null +++ b/src/refinement/refinement_pattern.jl @@ -0,0 +1,84 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int + + function CubicSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1 / sqrt(2), 1 / sqrt(2)] + direction_2 = [1 / sqrt(2), -1 / sqrt(2)] + direction_3 = -direction_1 + direction_4 = -direction_2 + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 4 + center_particle) + end +end + +struct TriangularSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int + + function TriangularSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + direction_1 = [0.0, 1.0] + direction_2 = [-sqrt(3) / 2, -0.5] + direction_3 = [sqrt(3) / 2, -0.5] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 3 + center_particle) + end +end + +struct HexagonalSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int + + function HexagonalSplitting(; epsilon=0.4, alpha=0.9, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1.0, 0.0] + direction_2 = [-1.0, 0.0] + direction_3 = [0.5, sqrt(3) / 2] + direction_4 = [0.5, -sqrt(3) / 2] + direction_5 = [-0.5, sqrt(3) / 2] + direction_6 = [-0.5, -sqrt(3) / 2] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4), + SVector{2, ELTYPE}(epsilon * direction_5), + SVector{2, ELTYPE}(epsilon * direction_6)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 6 + center_particle) + end +end diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl new file mode 100644 index 000000000..224517282 --- /dev/null +++ b/src/refinement/resize.jl @@ -0,0 +1,182 @@ +function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Resize all systems + foreach_system(semi) do system + capacity(system) > nparticles(system) && resize!(system, capacity(system)) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function Base.deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Delete at specific indices + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + deleteat!(system, v, u) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) + (; systems) = semi + + copyto!(_v_ode, v_ode) + copyto!(_u_ode, u_ode) + + # Get ranges after resizing the systems + ranges_v_new, ranges_u_new = ranges_vu(systems) + + ranges_v_old = copy(semi.ranges_v) + ranges_u_old = copy(semi.ranges_u) + + # Set ranges after resizing the systems + for i in 1:length(systems) + semi.ranges_v[i] = ranges_v_new[i] + semi.ranges_u[i] = ranges_u_new[i] + end + + for i in eachindex(ranges_u_old) + length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i])) + for j in 0:(length_u - 1) + u_ode[ranges_u_new[i][1] + j] = _u_ode[ranges_u_old[i][1] + j] + end + + length_v = min(length(ranges_v_old[i]), length(ranges_v_new[i])) + for j in 0:(length_v - 1) + v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j] + end + end + + sizes_u = sum(u_nvariables(system) * n_moving_particles(system) for system in systems) + sizes_v = sum(v_nvariables(system) * n_moving_particles(system) for system in systems) + + resize!(v_ode, sizes_v) + resize!(u_ode, sizes_u) + + resize!(_v_ode, sizes_v) + resize!(_u_ode, sizes_u) + + # TODO: Do the following in the callback + # resize!(integrator, (length(v_ode), length(u_ode))) + + # # Tell OrdinaryDiffEq that u has been modified + # u_modified!(integrator, true) + return v_ode +end + +Base.resize!(system::System, capacity_system) = system + +function Base.resize!(system::FluidSystem, capacity_system) + return resize!(system, system.particle_refinement, capacity_system) +end + +Base.resize!(system, ::Nothing, capacity_system) = system + +function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) + (; mass, pressure, density_calculator) = system + + resize!(mass, capacity_system) + resize!(pressure, capacity_system) + + resize_density!(system, capacity_system, density_calculator) + + resize_cache!(system, capacity_system) + + resize_refinement!(system) + + return system +end + +function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) + (; mass, density_calculator) = system + + resize!(mass, capacity_system) + + resize_density!(system, capacity_system, density_calculator) + + resize_cache!(system, capacity_system) + + resize_refinement!(system) + + return system +end + +resize_density!(system, n::Int, ::SummationDensity) = resize!(system.cache.density, n) +resize_density!(system, n::Int, ::ContinuityDensity) = system + +function resize_cache!(system::WeaklyCompressibleSPHSystem, n::Int) + resize!(system.cache.smoothing_length, n) + + return system +end + +function resize_cache!(system::EntropicallyDampedSPHSystem, n::Int) + resize!(system.cache.smoothing_length, n) + resize!(system.cache.beta, n) + resize!(system.cache.pressure_average, n) + resize!(system.cache.neighbor_counter, n) + + return system +end + +Base.deleteat!(system::System, v, u) = system + +function Base.deleteat!(system::FluidSystem, v, u) + return deleteat!(system, system.particle_refinement, v, u) +end + +Base.deleteat!(system::FluidSystem, ::Nothing, v, u) = system + +function Base.deleteat!(system::FluidSystem, refinement, v, u) + (; delete_candidates) = refinement + + delete_counter = 0 + + for particle in eachparticle(system) + if delete_candidates[particle] + # swap particles (keep -> delete) + dump_id = nparticles(system) - delete_counter + + vel_keep = current_velocity(v, system, dump_id) + pos_keep = current_coords(u, system, dump_id) + + mass_keep = hydrodynamic_mass(system, dump_id) + density_keep = particle_density(v, system, dump_id) + pressure_keep = particle_pressure(v, system, dump_id) + + system.cache.smoothing_length[particle] = smoothing_length(system, dump_id) + + system.mass[particle] = mass_keep + + set_particle_pressure!(v, system, particle, pressure_keep) + + set_particle_density!(v, system, particle, density_keep) + + for dim in 1:ndims(system) + v[dim, particle] = vel_keep[dim] + u[dim, particle] = pos_keep[dim] + end + + delete_counter += 1 + end + end + + resize!(system, nparticles(system) - delete_counter) + + return system +end + +@inline capacity(system) = nparticles(system) + +@inline capacity(system::FluidSystem) = capacity(system, system.particle_refinement) + +@inline capacity(system, ::Nothing) = nparticles(system) + +@inline function capacity(system, particle_refinement) + return particle_refinement.n_new_particles[] + nparticles(system) +end diff --git a/src/refinement/split.jl b/src/refinement/split.jl new file mode 100644 index 000000000..b7743d7d9 --- /dev/null +++ b/src/refinement/split.jl @@ -0,0 +1,98 @@ +function split_particles!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + collect_split_candidates!(system, v, u) + end + + resize!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + split_particles!(system, v, u) + end + + return semi +end + +@inline collect_split_candidates!(system, v, u) = System + +@inline function collect_split_candidates!(system::FluidSystem, v, u) + return collect_split_candidates!(system, system.particle_refinement, v, u) +end + +@inline collect_split_candidates!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) + (; mass_ref, max_spacing_ratio, split_candidates, refinement_pattern) = particle_refinement + (; n_children, center_particle) = refinement_pattern + + resize!(split_candidates, 0) + + for particle in eachparticle(system) + m_a = hydrodynamic_mass(system, particle) + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a > m_max + push!(split_candidates, particle) + end + end + + particle_refinement.n_new_particles[] = length(split_candidates) * + (n_children - center_particle) + + return system +end + +@inline split_particles!(system, v, u) = System + +@inline function split_particles!(system::FluidSystem, v, u) + return split_particles!(system, system.particle_refinement, v, u) +end + +@inline split_particles!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function split_particles!(system::FluidSystem, particle_refinement, v, u) + (; smoothing_length) = system.cache + (; split_candidates, refinement_pattern, n_new_particles) = particle_refinement + (; alpha, relative_position, n_children, center_particle) = refinement_pattern + + child_id_global = nparticles(system) - n_new_particles[] + for particle in split_candidates + smoothing_length_old = smoothing_length[particle] + mass_old = system.mass[particle] + + system.mass[particle] = mass_old / n_children + + p_a = particle_pressure(v, system, particle) + rho_a = particle_density(v, system, particle) + + smoothing_length[particle] = alpha * smoothing_length_old + + pos_center = current_coords(u, system, particle) + vel_center = current_velocity(v, system, particle) + + for child_id_local in 1:(n_children - center_particle) + child_id_global += 1 + + system.mass[child_id_global] = mass_old / n_children + + set_particle_pressure!(v, system, child_id_global, p_a) + + set_particle_density!(v, system, child_id_global, rho_a) + + smoothing_length[child_id_global] = alpha * smoothing_length_old + + rel_pos = smoothing_length_old * relative_position[child_id_local] + new_pos = pos_center + rel_pos + + for dim in 1:ndims(system) + u[dim, child_id_global] = new_pos[dim] + v[dim, child_id_global] = vel_center[dim] + end + end + end + + return system +end diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 8b7e507ed..dcf3d291a 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -69,6 +69,14 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, smoothing_length, viscosity, correction, cache) end +function smoothing_length(system, boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + +function smoothing_length(boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + @doc raw""" AdamiPressureExtrapolation(; pressure_offset=0.0) @@ -480,7 +488,7 @@ end resulting_acceleration = neighbor_system.acceleration - current_acceleration(system, particle) - kernel_weight = smoothing_kernel(boundary_model, distance) + kernel_weight = smoothing_kernel(boundary_model, distance, particle) pressure[particle] += (pressure_offset + @@ -581,12 +589,13 @@ end return density end -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) - (; smoothing_kernel, smoothing_length, correction) = system.boundary_model +# TODO +# @inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) +# (; smoothing_kernel, smoothing_length, correction) = system.boundary_model - return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, - smoothing_length, correction, system, particle) -end +# return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, +# smoothing_length, correction, system, particle) +# end @inline function correction_matrix(system::BoundarySystem, particle) extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 59f634092..d90d2a291 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -29,7 +29,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance, particle) continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff, grad_kernel, particle) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index a3f47f611..d37fed82e 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -328,6 +328,10 @@ end return kernel(smoothing_kernel, distance, smoothing_length) end +@inline function smoothing_length(system::BoundarySPHSystem, particle) + return smoothing_length(system, system.boundary_model, particle) +end + function update_positions!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t) (; movement) = system diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3644c3630..7e81ded7a 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,9 +30,10 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - dv_pressure = pressure_acceleration(particle_system, neighbor_system, neighbor, + dv_pressure = pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) @@ -44,20 +45,29 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term when using `TransportVelocityAdami` dv_convection = momentum_convection(particle_system, neighbor_system, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) + dv_velocity_correction = velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + + dv_velocity_correction[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - pressure_evolution!(dv, particle_system, v_diff, grad_kernel, - particle, pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) + pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, v_diff, grad_kernel, + particle, neighbor, pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) transport_velocity!(dv, particle_system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) @@ -69,11 +79,47 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end -@inline function pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, - pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) - (; smoothing_length) = particle_system +@inline function pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + v_diff, grad_kernel, particle, neighbor, + pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) + (; particle_refinement) = particle_system + + # This is basically the continuity equation times `sound_speed^2` + artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) + + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + dv[end, particle] += beta_inv_a * artificial_eos + + pressure_damping_term(particle_system, neighbor_system, + particle_refinement, particle, neighbor, + pos_diff, distance, beta_inv_a, m_a, m_b, + p_a, p_b, rho_b, rho_b, grad_kernel) + + pressure_reduction(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_inv_a, beta_inv_b) + + return dv +end +@inline beta_correction(system, particle) = one(eltype(system)) + +@inline function beta_correction(system::FluidSystem, particle) + beta_correction(system, system.particle_refinement, particle) +end + +@inline beta_correction(particle_system, ::Nothing, particle) = one(eltype(particle_system)) + +@inline function beta_correction(particle_system, refinement, particle) + return inv(particle_system.cache.beta[particle]) +end +function pressure_damping_term(particle_system, neighbor_system, ::Nothing, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -81,15 +127,13 @@ end # EDAC pressure evolution pressure_diff = p_a - p_b - # This is basically the continuity equation times `sound_speed^2` - artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) - eta_a = rho_a * particle_system.nu_edac eta_b = rho_b * particle_system.nu_edac eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) - # TODO For variable smoothing length use average smoothing length - tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length(particle_system, particle) + + smoothing_length(particle_system, particle)) + tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length_average^2) # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 # They argued that the formulation is more flexible because of the possibility to formulate @@ -102,11 +146,110 @@ end # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 # # This is similar to density diffusion in WCSPH - damping_term = volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) + return volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) +end - dv[end, particle] += artificial_eos + damping_term +function pressure_damping_term(particle_system, neighbor_system, refinement, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel_a) + (; sound_speed) = particle_system - return dv + # EDAC pressure evolution + pressure_diff = p_a - p_b + + # Haftu et al. (2022) uses `alpha_edac = 1.5` in all their simulations + alpha_edac = 1.5 + + # TODO: Haftu et al. (2022) use `8` but I think it depeneds on the dimension (see Monaghan, 2005) + tmp = 2 * ndims(particle_system) + 4 + + nu_edac_a = alpha_edac * sound_speed * smoothing_length(particle_system, particle) / tmp + nu_edac_b = alpha_edac * sound_speed * smoothing_length(neighbor_system, neighbor) / tmp + + nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) + + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + return beta_inv_a * nu_edac_ab * pressure_diff * dot(pos_diff, grad_W_avg) * m_b / rho_b +end + +function pressure_reduction(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + return zero(eltype(particle_system)) +end + +function pressure_reduction(particle_system, neighbor_system, refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = (p_a - p_a_avg) / (rho_a^2 * beta_inv_a) + P_b = (p_b - p_b_avg) / (rho_b^2 * beta_inv_b) + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + v_diff = advection_velocity(v_particle_system, particle_system, particle) - + current_velocity(v_particle_system, particle_system, particle) + + return m_b * (dot(v_diff, P_a * grad_kernel_a + P_b * grad_kernel_b)) +end + +@inline function velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + return zero(pos_diff) +end + +@inline function velocity_correction(particle_system, + neighbor_system::EntropicallyDampedSPHSystem, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + velocity_correction(particle_system, neighbor_system, + particle_system.particle_refinement, + pos_diff, distance, v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) +end + +@inline function velocity_correction(particle_system, neighbor_system, ::Nothing, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + return zero(pos_diff) +end + +@inline function velocity_correction(particle_system, neighbor_system, + particle_refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, particle_system, particle) + advection_velocity_a = advection_velocity(v_particle_system, particle_system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + v_diff = momentum_velocity_a - momentum_velocity_b + v_diff_tilde = advection_velocity_a - advection_velocity_b + + beta_inv_a = beta_correction(particle_system, particle) + + return -m_b * beta_inv_a * + (dot(v_diff_tilde - v_diff, grad_kernel) * momentum_velocity_a) / rho_b end # We need a separate method for EDAC since the density is stored in `v[end-1,:]`. @@ -125,3 +268,39 @@ end grad_kernel) return dv end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + return pressure_acceleration(particle_system, particle_system.particle_refinement, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, + correction) +end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + ::Nothing, particle, neighbor_system, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + (; pressure_acceleration_formulation) = particle_system + return pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) +end + +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + particle_refinement, neighbor_system, particle, + neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, + pos_diff, distance, W_a, correction) + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = beta_correction(particle_system, particle) * (p_a - p_a_avg) / rho_a^2 + P_b = beta_correction(neighbor_system, neighbor) * (p_b - p_b_avg) / rho_b^2 + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (P_a * grad_kernel_a + P_b * grad_kernel_b) +end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -44,12 +44,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K - smoothing_length :: ELTYPE sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -59,55 +58,57 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, transport_velocity :: TV source_terms :: ST buffer :: B + particle_refinement :: PR cache :: C +end - function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed; - pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), - transport_velocity=nothing, - alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel)), - source_terms=nothing, buffer_size=nothing) - buffer = isnothing(buffer_size) ? nothing : - SystemBuffer(nparticles(initial_condition), buffer_size) +function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, sound_speed; + pressure_acceleration=inter_particle_averaged_pressure, + density_calculator=SummationDensity(), + transport_velocity=nothing, + alpha=0.5, viscosity=nothing, + acceleration=ntuple(_ -> 0.0, + ndims(smoothing_kernel)), + particle_refinement=nothing, + source_terms=nothing, buffer_size=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) - initial_condition = allocate_buffer(initial_condition, buffer) + initial_condition = allocate_buffer(initial_condition, buffer) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) - mass = copy(initial_condition.mass) + mass = copy(initial_condition.mass) - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end + if ndims(smoothing_kernel) != NDIMS + throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) + end - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - nothing) + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, + density_calculator, + NDIMS, ELTYPE, + nothing) - nu_edac = (alpha * smoothing_length * sound_speed) / 8 + nu_edac = (alpha * smoothing_length * sound_speed) / 8 - cache = create_cache_density(initial_condition, density_calculator) - cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = create_cache_density(initial_condition, density_calculator) + cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) - new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), - typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms), - typeof(buffer), - typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, transport_velocity, source_terms, - buffer, cache) - end + return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, + smoothing_kernel, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, + transport_velocity, source_terms, buffer, + particle_refinement, cache) end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fd0a9c3a5..4a87c8646 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,8 +13,33 @@ function create_cache_density(ic, ::ContinuityDensity) return (;) end +function create_cache_refinement(initial_condition, ::Nothing, smoothing_length) + return (; smoothing_length) +end + +function create_cache_refinement(initial_condition, refinement, smoothing_length) + smoothng_length_factor = smoothing_length / initial_condition.particle_spacing + + beta = Vector{eltype(initial_condition)}(undef, nparticles(initial_condition)) + + return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)), + smoothing_length_factor=smoothng_length_factor, beta=beta) +end + @propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] +function smoothing_length(system::FluidSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::FluidSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::FluidSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + function write_u0!(u0, system::FluidSystem) (; initial_condition) = system @@ -89,7 +114,7 @@ include("entropically_damped_sph/entropically_damped_sph.jl") add_velocity!(du, v, particle, system, system.transport_velocity) end -@inline function momentum_convection(system, neighbor_system, +@inline function momentum_convection(system, neighbor_system, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return zero(grad_kernel) @@ -98,9 +123,59 @@ end @inline function momentum_convection(system, neighbor_system::Union{EntropicallyDampedSPHSystem, WeaklyCompressibleSPHSystem}, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) momentum_convection(system, neighbor_system, system.transport_velocity, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) + system.particle_refinement, pos_diff, distance, v_particle_system, + v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) +end + +function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + # Check if `UpdateCallback` is used when simulating with TVF + update_final!(system, system.transport_velocity, + v, u, v_ode, u_ode, semi, t; update_from_callback) + + # Compute correction factor when using particle refinement + compute_beta_correction!(system, system.particle_refinement, v_ode, u_ode, semi) + + return system +end + +compute_beta_correction!(system, ::Nothing, v_ode, u_ode, semi) = system + +function compute_beta_correction!(system, refinement, v_ode, u_ode, semi) + (; beta) = system.cache + + set_zero!(beta) + + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + + system_coords = current_coordinates(u, system) + + foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, + pos_diff, distance + rho_a = particle_density(v, system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + W_deriv = kernel_deriv(system.smoothing_kernel, distance, + smoothing_length(system, particle)) + + beta[particle] -= m_b * distance * W_deriv * (rho_a * ndims(system)) + end + end + + return system end diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index 5b0b82e11..62105cf6c 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -107,7 +107,7 @@ function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing end # Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction) (; pressure_acceleration_formulation) = particle_system @@ -118,7 +118,7 @@ end end # Formulation using asymmetric gradient formulation for corrections depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction::Union{KernelCorrection, diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 66140775e..9b736795e 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -105,10 +105,8 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, system_coords, neighbor_system_coords, neighborhood_search) 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 = smoothing_kernel_grad(system, pos_diff, distance, - particle) + density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] * smoothing_length diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bc8fdf005..374c81146 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -76,13 +76,14 @@ end return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) end -@inline function momentum_convection(system, neighbor_system, ::Nothing, +@inline function momentum_convection(system, neighbor_system, ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) end @inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) volume_a = m_a / rho_a @@ -101,6 +102,28 @@ end return volume_term * (0.5 * (A_a + A_b)) * grad_kernel end +@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, rho_a, rho_b, + m_a, m_b, particle, neighbor, grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, system, particle) + advection_velocity_a = advection_velocity(v_particle_system, system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + factor_a = beta_correction(system, particle) / rho_a + factor_b = beta_correction(neighbor_system, neighbor) / rho_b + + A_a = factor_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' + A_b = factor_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + + grad_kernel_a = smoothing_kernel_grad(system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (A_a * grad_kernel_a + A_b * grad_kernel_b) +end + @inline transport_velocity!(dv, system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) = dv @inline function transport_velocity!(dv, system::FluidSystem, @@ -154,12 +177,6 @@ function update_callback_used!(system, transport_velocity) system.cache.update_callback_used[] = true end -function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) - update_final!(system, system.transport_velocity, - v, u, v_ode, u_ode, semi, t; update_from_callback) -end - function update_final!(system::FluidSystem, ::Nothing, v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1a0665a71..7cef990a5 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -123,13 +123,20 @@ end v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) + # TODO do we need the average smoothing length here? pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length, grad_kernel, nu_a, nu_b) + smoothing_length_particle, grad_kernel, nu_a, nu_b) return m_b * pi_ab end @@ -170,12 +177,12 @@ end # Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". # In: Reports on Progress in Physics (2005), pages 1703-1759. # [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01) -function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan) - (; smoothing_length) = system +function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan, + smoothing_length_particle) (; alpha) = viscosity sound_speed = system_sound_speed(system) - return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4) + return alpha * smoothing_length_particle * sound_speed / (2 * ndims(system) + 4) end @doc raw""" @@ -213,13 +220,55 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system + viscosity(particle_system, neighbor_system, particle_system.particle_refinement, + v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +end + +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + beta_inv_a = beta_correction(particle_system, particle_refinement, particle) + + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length(particle_system, particle)) + + v_a = viscous_velocity(v_particle_system, particle_system, particle) + v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) + v_diff = v_a - v_b + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + tmp = (distance^2 + 0.001 * smoothing_length(particle_system, particle)^2) + + return m_b * beta_inv_a * 4 * nu_a * dot(pos_diff, grad_W_avg) * v_diff / + (tmp * (rho_a + rho_b)) +end +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -230,8 +279,8 @@ end eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - # TODO For variable smoothing_length use average smoothing length - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length_particle + smoothing_length_neighbor) + tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) volume_a = m_a / rho_a volume_b = m_b / rho_b @@ -250,7 +299,7 @@ end return visc .* v_diff end -function kinematic_viscosity(system, viscosity::ViscosityAdami) +function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length) return viscosity.nu end diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c10740ab8..33a1b7787 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -75,9 +75,9 @@ end @inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b, pos_diff, distance, system, particle, neighbor) - (; smoothing_length) = system - - return ((rho_a - rho_b) / 2smoothing_length) * pos_diff / distance + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(system, particle) + return ((rho_a - rho_b) / (2 * smoothing_length_particle)) * pos_diff / distance end @doc raw""" @@ -207,7 +207,7 @@ end distance < sqrt(eps(typeof(distance))) && return (; delta) = density_diffusion - (; smoothing_length, state_equation) = particle_system + (; state_equation) = particle_system (; sound_speed) = state_equation volume_b = m_b / rho_b @@ -216,7 +216,10 @@ end particle_system, particle, neighbor) density_diffusion_term = dot(psi, grad_kernel) * volume_b - dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(particle_system, particle) + dv[end, particle] += delta * smoothing_length_particle * sound_speed * + density_diffusion_term end # Density diffusion `nothing` or interaction other than fluid-fluid diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0fe58ab21..b874737bf 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -53,7 +53,8 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor) dv_pressure = pressure_correction * - pressure_acceleration(particle_system, neighbor_system, neighbor, + pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 6d66ba06a..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,15 +44,14 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} +struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, + PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} density_calculator :: DC state_equation :: SE smoothing_kernel :: K - smoothing_length :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V density_diffusion :: DD @@ -62,6 +61,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, source_terms :: ST surface_tension :: SRFT buffer :: B + particle_refinement :: PR # TODO cache :: C end @@ -76,6 +76,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, + particle_refinement=nothing, surface_tension=nothing) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -116,14 +117,18 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, - source_terms, surface_tension, buffer, cache) + source_terms, surface_tension, buffer, + particle_refinement, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..114f72aed 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -32,7 +32,7 @@ end rho_b = neighbor_system.material_density[neighbor] grad_kernel = smoothing_kernel_grad(particle_system, initial_pos_diff, - initial_distance) + initial_distance, particle) m_a = particle_system.mass[particle] m_b = neighbor_system.mass[neighbor] @@ -88,7 +88,7 @@ function interact!(dv, v_particle_system, u_particle_system, # 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) + grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance, particle) # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles # corresponding to the chosen boundary model. @@ -99,7 +99,8 @@ function interact!(dv, v_particle_system, u_particle_system, # 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, + dv_boundary = pressure_acceleration(neighbor_system, particle_system, + neighbor, particle, m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, distance, grad_kernel, neighbor_system.correction) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index dfd71a2d6..8275da22b 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -286,7 +286,7 @@ end pos_diff = current_coords(system, particle) - current_coords(system, neighbor) grad_kernel = smoothing_kernel_grad(system, initial_pos_diff, - initial_distance) + initial_distance, particle) result = volume * pos_diff * grad_kernel' diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ef5f68cc6..90c5fa435 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -252,7 +252,7 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["viscosity"] = type2string(system.viscosity) write2vtk!(vtk, system.viscosity) vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length"] = system.smoothing_length + # vtk["smoothing_length"] = system.smoothing_length TODO vtk["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem