From 226e23ffbea40d8f55a90de4d9fd794e9c8fbd1f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 16:38:58 +0100 Subject: [PATCH 01/27] add resize --- src/TrixiParticles.jl | 1 + src/multi_resolution/multi_resolution.jl | 1 + src/multi_resolution/resize.jl | 170 +++++++++++++++++++++++ 3 files changed, 172 insertions(+) create mode 100644 src/multi_resolution/multi_resolution.jl create mode 100644 src/multi_resolution/resize.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bcdce98ac..21dbdd8a7 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -40,6 +40,7 @@ include("general/system.jl") # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded` include("util.jl") include("preprocessing/preprocessing.jl") +include("multi_resolution/multi_resolution.jl") include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl new file mode 100644 index 000000000..648a9fb1f --- /dev/null +++ b/src/multi_resolution/multi_resolution.jl @@ -0,0 +1 @@ +include("resize.jl") diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl new file mode 100644 index 000000000..23203e5c4 --- /dev/null +++ b/src/multi_resolution/resize.jl @@ -0,0 +1,170 @@ +function resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Resize all systems + foreach_system(semi) do system + resize!(system, capacity(system)) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function 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 resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) + copyto!(_v_ode, v_ode) + copyto!(_u_ode, u_ode) + + # Get ranges after resizing the systems + ranges_v_new, ranges_u_new = ranges_vu(semi.systems) + + ranges_v_old = semi.ranges_v + ranges_u_old = semi.ranges_u + + # Set ranges after resizing the systems + for i in 1:length(semi.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 1:length_u + 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 1:length_v + v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j] + end + end + + capacity_global = sum(system -> nparticles(system), semi.systems) + + resize!(v_ode, capacity_global) + resize!(u_ode, capacity_global) + + resize!(_v_ode, capacity_global) + resize!(_u_ode, capacity_global) + + # 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 + +resize!(system, capacity_system) = system + +function resize!(system::FluidSystem, capacity_system) + return resize!(system, system.particle_refinement, capacity_system) +end + +resize!(system, ::Nothing, capacity_system) = system + +function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) + (; mass, pressure, cache, density_calculator) = system + + refinement.n_particles_before_resize = nparticles(system) + + resize!(mass, capacity_system) + resize!(pressure, capacity_system) + resize_density!(system, capacity_system, density_calculator) + resize_cache!(system, cache, n) +end + +function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) + (; mass, cache, density_calculator) = system + + refinement.n_particles_before_resize = nparticles(system) + + resize!(mass, capacity_system) + resize_density!(system, capacity_system, density_calculator) + resize_cache!(system, capacity_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, n::Int) + resize!(system.cache.smoothing_length, n) + + return system +end + +function resize_cache!(system::EntropicallyDampedSPHSystem, n) + resize!(system.cache.smoothing_length, n) + resize!(system.cache.pressure_average, n) + resize!(system.cache.neighbor_counter, n) + + return system +end + +deleteat!(system, v, u) = system + +function deleteat!(system::FluidSystem, v, u) + return deleteat!(system, system.particle_refinement, v, u) +end + +deleteat!(system, ::Nothing, v, u) = system + +function deleteat!(system::FluidSystem, refinement, v, u) + (; delete_candidates) = refinement + + delete_counter = 0 + + for particle in eachparticle(system) + if !iszero(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(system, v, dump_id) + pressure_keep = particle_pressure(system, v, dump_id) + smoothing_length_keep = smoothing_length(system, dump_id) + + system.mass[particle] = mass_keep + system.cache.smoothing_length[particle] = smoothing_length_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) = 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 From afe21961f87772377937257d2b4cb9bc3e55fdc4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 17:15:23 +0100 Subject: [PATCH 02/27] adapt for tests --- src/TrixiParticles.jl | 2 +- src/multi_resolution/multi_resolution.jl | 1 + src/multi_resolution/particle_refinement.jl | 8 ++ src/multi_resolution/resize.jl | 12 +-- .../fluid/entropically_damped_sph/system.jl | 78 +++++++++---------- 5 files changed, 56 insertions(+), 45 deletions(-) create mode 100644 src/multi_resolution/particle_refinement.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 21dbdd8a7..aa685c70d 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -40,7 +40,6 @@ include("general/system.jl") # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded` include("util.jl") include("preprocessing/preprocessing.jl") -include("multi_resolution/multi_resolution.jl") include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") @@ -52,6 +51,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("multi_resolution/multi_resolution.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl index 648a9fb1f..547ee3c35 100644 --- a/src/multi_resolution/multi_resolution.jl +++ b/src/multi_resolution/multi_resolution.jl @@ -1 +1,2 @@ include("resize.jl") +include("particle_refinement.jl") diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl new file mode 100644 index 000000000..b5d6f8340 --- /dev/null +++ b/src/multi_resolution/particle_refinement.jl @@ -0,0 +1,8 @@ +struct ParticleRefinement + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} +end + +function ParticleRefinement() + return ParticleRefinement(Ref(0), Ref(0)) +end diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index 23203e5c4..a32556c30 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -77,22 +77,24 @@ resize!(system, ::Nothing, capacity_system) = system function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) (; mass, pressure, cache, density_calculator) = system - refinement.n_particles_before_resize = nparticles(system) + refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize!(pressure, capacity_system) resize_density!(system, capacity_system, density_calculator) - resize_cache!(system, cache, n) + # TODO + # resize_cache!(system, cache, n) end function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) (; mass, cache, density_calculator) = system - refinement.n_particles_before_resize = nparticles(system) + refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize_density!(system, capacity_system, density_calculator) - resize_cache!(system, capacity_system) + # TODO + # resize_cache!(system, capacity_system) return system end @@ -166,5 +168,5 @@ end @inline capacity(system, ::Nothing) = nparticles(system) @inline function capacity(system, particle_refinement) - return particle_refinement.n_new_particles + nparticles(system) + return particle_refinement.n_new_particles[] + nparticles(system) end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..94064555d 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -44,7 +44,7 @@ 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 @@ -59,55 +59,55 @@ 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) / + (2 * ndims(initial_condition) + 4) - 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...) - 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, smoothing_length, 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) From 4162c29d862cc0ed003d31658256f9928d228199 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 15 Nov 2024 14:44:25 +0100 Subject: [PATCH 03/27] fix `resize!(system, n)` --- src/general/semidiscretization.jl | 19 +++++++++++-------- src/multi_resolution/resize.jl | 26 ++++++++++++++------------ 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..cf9df1b5e 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 diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index a32556c30..2c7fa9405 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -1,4 +1,4 @@ -function resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) +function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) # Resize all systems foreach_system(semi) do system resize!(system, capacity(system)) @@ -22,15 +22,15 @@ function deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) return semi end -function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) +function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) copyto!(_v_ode, v_ode) copyto!(_u_ode, u_ode) # Get ranges after resizing the systems ranges_v_new, ranges_u_new = ranges_vu(semi.systems) - ranges_v_old = semi.ranges_v - ranges_u_old = semi.ranges_u + 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(semi.systems) @@ -40,12 +40,12 @@ function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) for i in eachindex(ranges_u_old) length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i])) - for j in 1:length_u + 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 1:length_v + 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 @@ -66,15 +66,15 @@ function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) return v_ode end -resize!(system, capacity_system) = system +Base.resize!(system::System, capacity_system) = system -function resize!(system::FluidSystem, capacity_system) +function Base.resize!(system::FluidSystem, capacity_system) return resize!(system, system.particle_refinement, capacity_system) end -resize!(system, ::Nothing, capacity_system) = system +Base.resize!(system, ::Nothing, capacity_system) = system -function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) +function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) (; mass, pressure, cache, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) @@ -86,7 +86,7 @@ function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_syste # resize_cache!(system, cache, n) end -function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) +function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) (; mass, cache, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) @@ -163,7 +163,9 @@ function deleteat!(system::FluidSystem, refinement, v, u) return system end -@inline capacity(system) = capacity(system, system.particle_refinement) +@inline capacity(system) = nparticles(system) + +@inline capacity(system::FluidSystem) = capacity(system, system.particle_refinement) @inline capacity(system, ::Nothing) = nparticles(system) From 9c0d0beb0d801ab5083ca6354366ffdda773cf6b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 15 Nov 2024 14:56:15 +0100 Subject: [PATCH 04/27] fix `deleteat!` --- src/multi_resolution/particle_refinement.jl | 3 ++- src/multi_resolution/resize.jl | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl index b5d6f8340..bf886ac12 100644 --- a/src/multi_resolution/particle_refinement.jl +++ b/src/multi_resolution/particle_refinement.jl @@ -1,8 +1,9 @@ struct ParticleRefinement n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles end function ParticleRefinement() - return ParticleRefinement(Ref(0), Ref(0)) + return ParticleRefinement(Ref(0), Ref(0), Bool[]) end diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index 2c7fa9405..87dc97fa5 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -9,7 +9,7 @@ function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) return semi end -function deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) +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) @@ -116,15 +116,15 @@ function resize_cache!(system::EntropicallyDampedSPHSystem, n) return system end -deleteat!(system, v, u) = system +Base.deleteat!(system::System, v, u) = system -function deleteat!(system::FluidSystem, v, u) +function Base.deleteat!(system::FluidSystem, v, u) return deleteat!(system, system.particle_refinement, v, u) end -deleteat!(system, ::Nothing, v, u) = system +Base.deleteat!(system::FluidSystem, ::Nothing, v, u) = system -function deleteat!(system::FluidSystem, refinement, v, u) +function Base.deleteat!(system::FluidSystem, refinement, v, u) (; delete_candidates) = refinement delete_counter = 0 @@ -138,12 +138,13 @@ function deleteat!(system::FluidSystem, refinement, v, u) pos_keep = current_coords(u, system, dump_id) mass_keep = hydrodynamic_mass(system, dump_id) - density_keep = particle_density(system, v, dump_id) - pressure_keep = particle_pressure(system, v, dump_id) - smoothing_length_keep = smoothing_length(system, dump_id) + density_keep = particle_density(v, system, dump_id) + pressure_keep = particle_pressure(v, system,dump_id) + #TODO + # smoothing_length_keep = smoothing_length(system, dump_id) + # system.cache.smoothing_length[particle] = smoothing_length_keep system.mass[particle] = mass_keep - system.cache.smoothing_length[particle] = smoothing_length_keep set_particle_pressure!(v, system, particle, pressure_keep) From 520835cb256657ef79964181272f0950347d2c69 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:24:47 +0100 Subject: [PATCH 05/27] restructure --- src/TrixiParticles.jl | 2 +- src/multi_resolution/multi_resolution.jl | 2 -- .../particle_refinement.jl => refinement/refinement.jl} | 2 ++ src/{multi_resolution => refinement}/resize.jl | 0 4 files changed, 3 insertions(+), 3 deletions(-) delete mode 100644 src/multi_resolution/multi_resolution.jl rename src/{multi_resolution/particle_refinement.jl => refinement/refinement.jl} (93%) rename src/{multi_resolution => refinement}/resize.jl (100%) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index aa685c70d..30594da02 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -51,7 +51,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") -include("multi_resolution/multi_resolution.jl") +include("refinement/refinement.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl deleted file mode 100644 index 547ee3c35..000000000 --- a/src/multi_resolution/multi_resolution.jl +++ /dev/null @@ -1,2 +0,0 @@ -include("resize.jl") -include("particle_refinement.jl") diff --git a/src/multi_resolution/particle_refinement.jl b/src/refinement/refinement.jl similarity index 93% rename from src/multi_resolution/particle_refinement.jl rename to src/refinement/refinement.jl index bf886ac12..87f06afef 100644 --- a/src/multi_resolution/particle_refinement.jl +++ b/src/refinement/refinement.jl @@ -1,3 +1,5 @@ +include("resize.jl") + struct ParticleRefinement n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} diff --git a/src/multi_resolution/resize.jl b/src/refinement/resize.jl similarity index 100% rename from src/multi_resolution/resize.jl rename to src/refinement/resize.jl From ea11f8438733791e88148f3c77c9347d3745c188 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:33:43 +0100 Subject: [PATCH 06/27] include order --- src/TrixiParticles.jl | 2 +- src/refinement/refinement.jl | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 30594da02..89026040c 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -51,7 +51,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") -include("refinement/refinement.jl") +include("refinement/resize.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 87f06afef..bf886ac12 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,5 +1,3 @@ -include("resize.jl") - struct ParticleRefinement n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} From 25cfff804da4d0153c82a610f2860f2070fcb17c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:38:40 +0100 Subject: [PATCH 07/27] apply formatter --- src/refinement/resize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index 87dc97fa5..a7e70fa47 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -139,7 +139,7 @@ function Base.deleteat!(system::FluidSystem, refinement, v, u) mass_keep = hydrodynamic_mass(system, dump_id) density_keep = particle_density(v, system, dump_id) - pressure_keep = particle_pressure(v, system,dump_id) + pressure_keep = particle_pressure(v, system, dump_id) #TODO # smoothing_length_keep = smoothing_length(system, dump_id) # system.cache.smoothing_length[particle] = smoothing_length_keep From 679d297481203fbe56e83c625d48c6974cea6d71 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 14 Nov 2024 16:21:41 +0100 Subject: [PATCH 08/27] Boundary Pressure Calculation for FSI (#498) * change how adami pressure extrapolation is calculated and add optional offset * remove unused function * update * move to dispatch on function * fix * fix * format * fix test * format * change how adami pressure extrapolation is calculated and add optional offset * remove unused function * update * move to dispatch on function * fix * fix * format * fix test * format * fix merge * fix merge * fix * fix * fix * add new bnd density calculator * missing code * fix test * fix * fix * fix docs * fix * format * review comments * fix test * fix test * format * fix * update docs * fix * set test up for 1.11 * format * implement suggestions * fix equation * formatting * format * Increase errors for 1.11 * Fix invalidations * Fix tests * Update ci.yml * revert * Update ci.yml * Update test/validation/validation.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * review comments * adjust docs * adjust docs * format * add complete equation * Update docs/src/systems/boundary.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- docs/src/systems/boundary.md | 31 ++- examples/fsi/falling_spheres_2d.jl | 2 +- src/TrixiParticles.jl | 4 +- .../dummy_particles/dummy_particles.jl | 181 +++++++++++++----- .../dummy_particles/dummy_particles.jl | 24 ++- 5 files changed, 175 insertions(+), 67 deletions(-) diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 29384b6d1..7beb184ed 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -55,13 +55,17 @@ of the boundary particle ``b``. ### Hydrodynamic density of dummy particles -We provide five options to compute the boundary density and pressure, determined by the `density_calculator`: +We provide six options to compute the boundary density and pressure, determined by the `density_calculator`: 1. (Recommended) With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the fluid according to [Adami et al., 2012](@cite Adami2012), and the density is obtained by applying the inverse of the state equation. This option usually yields the best results of the options listed here. -2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, +2. (Only relevant for FSI) With [`BernoulliPressureExtrapolation`](@ref), the pressure is extrapolated from the + pressure similar to the [`AdamiPressureExtrapolation`](@ref), but a relative velocity-dependent pressure part + is calculated between moving solids and fluids, which increases the boundary pressure in areas prone to + penetrations. +3. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, and the pressure is computed from the density with the state equation. -3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, +4. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, and the pressure is computed from the density with the state equation. Note that this causes a gap between fluid and boundary where the boundary is initialized without any contact to the fluid. This is due to overestimation of the boundary density @@ -69,10 +73,10 @@ We provide five options to compute the boundary density and pressure, determined contact to the fluid. Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat. See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13). -4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure +5. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure is computed from the density with the state equation. This option is not recommended. The other options yield significantly better results. -5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure +6. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure is not used. Instead, the fluid pressure is mirrored as boundary pressure in the momentum equation. This option is not recommended due to stability issues. See [`PressureMirroring`](@ref) @@ -93,7 +97,20 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den AdamiPressureExtrapolation ``` -#### 4. [`PressureZeroing`](@ref) +#### 2. [`BernoulliPressureExtrapolation`](@ref) +Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation: +```math +p_b = \frac{\sum_f (p_f + \frac{1}{2} \, \rho_{\text{neighbor}} \left( \frac{ (\mathbf{v}_f - \mathbf{v}_{\text{body}}) \cdot (\mathbf{x}_f - \mathbf{x}_{\text{neighbor}}) }{ \left\| \mathbf{x}_f - \mathbf{x}_{\text{neighbor}} \right\| } \right)^2 \times \text{factor} +\rho_f (\bm{g} - \bm{a}_b) \cdot \bm{r}_{bf}) W(\Vert r_{bf} \Vert, h)}{\sum_f W(\Vert r_{bf} \Vert, h)} +``` +where ``\mathbf{v}_f`` is the velocity of the fluid and ``\mathbf{v}_{\text{body}}`` is the velocity of the body. +This adjustment provides a higher boundary pressure for solid bodies moving with a relative velocity to the fluid to prevent penetration. +This modification is original and not derived from any literature source. + +```@docs + BernoulliPressureExtrapolation +``` + +#### 5. [`PressureZeroing`](@ref) This is the simplest way to implement dummy boundary particles. The density of each particle is set to the reference density and the pressure to the @@ -102,7 +119,7 @@ reference pressure (the corresponding pressure to the reference density by the s PressureZeroing ``` -#### 5. [`PressureMirroring`](@ref) +#### 6. [`PressureMirroring`](@ref) Instead of calculating density and pressure for each boundary particle, we modify the momentum equation, diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 6a2397a4b..a06056f8f 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -63,7 +63,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary -boundary_density_calculator = AdamiPressureExtrapolation() +boundary_density_calculator = BernoulliPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 89026040c..7f2d2272b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -70,7 +70,9 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka + PressureMirroring, PressureZeroing, BoundaryModelLastiwka, + BernoulliPressureExtrapolation + export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 1cd8e0dfe..8b7e507ed 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -70,11 +70,43 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, end @doc raw""" - AdamiPressureExtrapolation() + AdamiPressureExtrapolation(; pressure_offset=0.0) `density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. + """ -struct AdamiPressureExtrapolation end +struct AdamiPressureExtrapolation{ELTYPE} + pressure_offset::ELTYPE + + function AdamiPressureExtrapolation(; pressure_offset=0.0) + return new{eltype(pressure_offset)}(pressure_offset) + end +end + +@doc raw""" + BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=1.0) + +`density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. +- `factor=1.0` : Setting `factor` allows to just increase the strength of the dynamic + pressure part. + +""" +struct BernoulliPressureExtrapolation{ELTYPE} + pressure_offset :: ELTYPE + factor :: ELTYPE + + function BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=0.0) + return new{eltype(pressure_offset)}(pressure_offset, factor) + end +end @doc raw""" PressureMirroring() @@ -134,7 +166,9 @@ end @inline create_cache_model(initial_density, ::ContinuityDensity) = (; initial_density) -function create_cache_model(initial_density, ::AdamiPressureExtrapolation) +function create_cache_model(initial_density, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}) density = copy(initial_density) volume = similar(initial_density) @@ -194,7 +228,7 @@ end # Note that the other density calculators are dispatched in `density_calculators.jl` @inline function particle_density(v, ::Union{AdamiPressureExtrapolation, PressureMirroring, - PressureZeroing}, + PressureZeroing, BernoulliPressureExtrapolation}, boundary_model, particle) (; cache) = boundary_model @@ -216,10 +250,11 @@ end function compute_density!(boundary_model, ::Union{ContinuityDensity, AdamiPressureExtrapolation, + BernoulliPressureExtrapolation, PressureMirroring, PressureZeroing}, system, v, u, v_ode, u_ode, semi) # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`. - # For `AdamiPressureExtrapolation`, the density is updated in `compute_pressure!`. + # For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. return boundary_model end @@ -299,8 +334,10 @@ end boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) end -function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, - v_ode, u_ode, semi) +function compute_pressure!(boundary_model, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}, + system, v, u, v_ode, u_ode, semi) (; pressure, cache, viscosity) = boundary_model set_zero!(pressure) @@ -327,16 +364,18 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, nhs = get_neighborhood_search(neighbor_system, system, semi) # Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) else nhs = get_neighborhood_search(system, neighbor_system, semi) # Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) end @threaded system for particle in eachparticle(system) @@ -378,67 +417,81 @@ function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZe return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system, system_coords, - neighbor_coords, v_neighbor_system, - neighborhood_search) +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, system_coords, + neighbor_coords, v, + v_neighbor_system, + neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, - neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, + v, v_neighbor_system, + neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator foreach_point_neighbor(neighbor_system, system, neighbor_coords, system_coords, neighborhood_search; points=eachparticle(neighbor_system), parallel=false) do neighbor, particle, pos_diff, distance # Since neighbor and particle are switched pos_diff = -pos_diff - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) +@inline function boundary_pressure_extrapolation!(boundary_model, system, neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator # Loop over all pairs of particles and neighbors within the kernel cutoff. foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, neighborhood_search; points=eachparticle(system)) do particle, neighbor, pos_diff, distance - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) +@inline function boundary_pressure_inner!(boundary_model, boundary_density_calculator, + system, neighbor_system::FluidSystem, v, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure, + pressure_offset) density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) - resulting_acc = neighbor_system.acceleration - - current_acceleration(system, particle) + resulting_acceleration = neighbor_system.acceleration - + current_acceleration(system, particle) kernel_weight = smoothing_kernel(boundary_model, distance) - pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, - neighbor) + - dot(resulting_acc, density_neighbor * pos_diff)) * + pressure[particle] += (pressure_offset + + + particle_pressure(v_neighbor_system, neighbor_system, + neighbor) + + + dynamic_pressure(boundary_density_calculator, density_neighbor, + v, v_neighbor_system, pos_diff, distance, + particle, neighbor, system, neighbor_system) + + + dot(resulting_acceleration, density_neighbor * pos_diff)) * kernel_weight cache.volume[particle] += kernel_weight @@ -447,14 +500,46 @@ end kernel_weight, particle, neighbor) end +@inline function dynamic_pressure(boundary_density_calculator, density_neighbor, v, + v_neighbor_system, pos_diff, distance, particle, neighbor, + system, neighbor_system) + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::BoundarySystem, neighbor_system) + if system.ismoving[] + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) + + return 0.5 * boundary_density_calculator.factor * density_neighbor * + normal_velocity^2 / distance + end + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::SolidSystem, neighbor_system) + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) / distance + + return boundary_density_calculator.factor * density_neighbor * + dot(normal_velocity, normal_velocity) / 2 +end + function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) return cache end -function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, - neighbor_system, v_neighbor_system, kernel_weight, - particle, neighbor) +function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, neighbor_system, + v_neighbor_system, kernel_weight, particle, neighbor) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) for dim in 1:ndims(neighbor_system) diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 27caa0293..a72ac409c 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -54,12 +54,14 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, - v_fluid .* - ones(size(fluid.coordinates)), - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid .* + ones(size(fluid.coordinates)), + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps() @@ -92,10 +94,12 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, v_fluid, - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid, + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps() From fea431cecad3d1de010fbbbd94202e5b5b3e9433 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 14 Nov 2024 18:41:58 +0100 Subject: [PATCH 09/27] Test docs on win (#664) * Test docs on win * Fix error on windows * Format * Update Documenter.yml * Update Documenter.yml --- .github/workflows/Documenter.yml | 9 ++++++++- docs/make.jl | 5 ++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index a748be851..5e093c7b7 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -25,7 +25,10 @@ concurrency: jobs: build-docs: name: Build and Deploy Documentation - runs-on: ubuntu-latest + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, windows-latest] steps: - name: Check out project uses: actions/checkout@v4 @@ -40,7 +43,11 @@ jobs: - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy + if: matrix.os == 'ubuntu-latest' env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key run: julia --project=docs --color=yes docs/make.jl + - name: Just Build + if: matrix.os != 'ubuntu-latest' + run: julia --project=docs --color=yes docs/make.jl diff --git a/docs/make.jl b/docs/make.jl index 3e3411f14..82d76064b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,9 +19,12 @@ function copy_file(filename, replaces...; content = read(source_path, String) content = replace(content, replaces...) + # Use `replace` to make sure the path uses forward slashes for URLs + filename_url = replace(filename, "\\" => "/") + header = """ ```@meta - EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename" + EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename_url" ``` """ content = header * content From fc3a38d78a146755d3532b89ec9d741deeb26804 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 14 Nov 2024 19:37:48 +0100 Subject: [PATCH 10/27] Add docs for GPU support (#660) * Add docs for GPU support * Reformat code * Implement suggestions * Implement suggestions * Fix typo --- docs/make.jl | 3 +- docs/src/gpu.md | 61 ++++++++++++++++++++++++++++ docs/src/overview.md | 8 +++- docs/src/reference-pointneighbors.md | 2 +- test/examples/examples.jl | 12 ++++++ 5 files changed, 83 insertions(+), 3 deletions(-) create mode 100644 docs/src/gpu.md diff --git a/docs/make.jl b/docs/make.jl index 82d76064b..ac541b0a1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -120,7 +120,8 @@ makedocs(sitename="TrixiParticles.jl", "Preprocessing" => [ "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md") ], - "Components" => [ + "GPU Support" => "gpu.md", + "API Reference" => [ "Overview" => "overview.md", "General" => [ "Semidiscretization" => joinpath("general", "semidiscretization.md"), diff --git a/docs/src/gpu.md b/docs/src/gpu.md new file mode 100644 index 000000000..8e8a12708 --- /dev/null +++ b/docs/src/gpu.md @@ -0,0 +1,61 @@ +# GPU Support + +GPU support is still an experimental feature that is actively being worked on. +As of now, the [`WeaklyCompressibleSPHSystem`](@ref) and the [`BoundarySPHSystem`](@ref) +are supported on GPUs. +We have tested this on GPUs by Nvidia and AMD. + +To run a simulation on a GPU, we need to use the [`FullGridCellList`](@ref) +as cell list for the [`GridNeighborhoodSearch`](@ref). +This cell list requires a bounding box for the domain, unlike the default cell list, which +uses an unbounded domain. +For simulations that are bounded by a closed tank, we can use the boundary of the tank +to obtain the bounding box as follows. +```jldoctest gpu; output=false, setup=:(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing)) +search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length) +min_corner = minimum(tank.boundary.coordinates, dims=2) .- search_radius +max_corner = maximum(tank.boundary.coordinates, dims=2) .+ search_radius +cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, max_corner) + +# output +PointNeighbors.FullGridCellList{PointNeighbors.DynamicVectorOfVectors{Int32, Matrix{Int32}, Vector{Int32}, Base.RefValue{Int32}}, Nothing, SVector{2, Float64}, SVector{2, Float64}}(Vector{Int32}[], nothing, [-0.24500000000000002, -0.24500000000000002], [1.245, 1.245]) +``` + +We then need to pass this cell list to the neighborhood search and the neighborhood search +to the [`Semidiscretization`](@ref). +```jldoctest gpu; output=false +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(; cell_list)) + +# output +┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ +│ Semidiscretization │ +│ ══════════════════ │ +│ #spatial dimensions: ………………………… 2 │ +│ #systems: ……………………………………………………… 2 │ +│ neighborhood search: ………………………… GridNeighborhoodSearch │ +│ total #particles: ………………………………… 636 │ +└──────────────────────────────────────────────────────────────────────────────────────────────────┘ +``` + +At this point, we should run the simulation and make sure that it still works and that +the bounding box is large enough. +For some simulations where particles move outside the initial tank coordinates, +for example when the tank is not closed or when the tank is moving, an appropriate +bounding box has to be specified. + +Then, we only need to specify the data type that is used for the simulation. +On an Nvidia GPU, we specify: +```julia +using CUDA +ode = semidiscretize(semi, tspan, data_type=CuArray) +``` +On an AMD GPU, we use: +```julia +using AMDGPU +ode = semidiscretize(semi, tspan, data_type=ROCArray) +``` +Then, we can run the simulation as usual. +All data is transferred to the GPU during initialization and all loops over particles +and their neighbors will be executed on the GPU as kernels generated by KernelAbstractions.jl. +Data is only copied to the CPU for saving VTK files via the [`SolutionSavingCallback`](@ref). diff --git a/docs/src/overview.md b/docs/src/overview.md index a113d1f9f..5d93c3732 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -1,10 +1,16 @@ # Overview +The actual API reference is not listed on a single page, like in most Julia packages, +but instead is split into multiple sections that follow a similar structure +as the code files themselves. +In these sections, API docs are combined with explanations of the theoretical background +of these methods. + The following page gives a rough overview of important parts of the code. ## Program flow To initiate a simulation, the goal is to solve an ordinary differential equation, for example, -by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then +by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then utilized to integrate ``\mathrm{d}u/\mathrm{d}t`` and ``\mathrm{d}v/\mathrm{d}t``, where ``u`` represents the particles' positions and ``v`` their properties such as velocity and density. During a single time step or an intermediate step of the time integration scheme, the functions diff --git a/docs/src/reference-pointneighbors.md b/docs/src/reference-pointneighbors.md index f8f8c5ac7..ae184bb2a 100644 --- a/docs/src/reference-pointneighbors.md +++ b/docs/src/reference-pointneighbors.md @@ -1,4 +1,4 @@ -# PointNeighbors.jl API +# [PointNeighbors.jl API](@id pointneighbors) ```@meta CurrentModule = PointNeighbors diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 5f42aa230..2e5462a35 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -9,8 +9,20 @@ "hydrostatic_water_column_2d.jl"), fluid_system=nothing, sol=nothing, semi=nothing, ode=nothing) + # Neighborhood search for `FullGridCellList` test below + search_radius = TrixiParticles.compact_support(smoothing_kernel, + smoothing_length) + min_corner = minimum(tank.boundary.coordinates, dims=2) .- search_radius + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ search_radius + cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, + max_corner) + semi_fullgrid = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(; + cell_list)) + hydrostatic_water_column_tests = Dict( "WCSPH default" => (), + "WCSPH with FullGridCellList" => (semi=semi_fullgrid,), "WCSPH with source term damping" => (source_terms=SourceTermDamping(damping_coefficient=1e-4),), "WCSPH with SummationDensity" => (fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), From 230aafce61f29d1a9cdfb5bc07e31251e019a84a Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 15 Nov 2024 15:49:46 +0100 Subject: [PATCH 11/27] Avoid bounds checking where it is safe to do so (#656) * Avoid bounds checking where it is safe to do so * Avoid one more bounds check in density diffusion * Reformat code * Use `extract_svector` from PointNeighbors.jl * Revert 09ab7bac2f8e3c95fc37fdf6809503d5362383e0 * Fix tests --- src/TrixiParticles.jl | 1 + src/general/density_calculators.jl | 6 +- src/general/system.jl | 24 ++++--- src/schemes/fluid/fluid.jl | 2 +- src/schemes/fluid/viscosity.jl | 36 ++++++---- .../density_diffusion.jl | 8 +-- .../fluid/weakly_compressible_sph/rhs.jl | 70 +++++++++++-------- .../fluid/weakly_compressible_sph/system.jl | 3 +- test/systems/solid_system.jl | 2 +- 9 files changed, 87 insertions(+), 65 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 7f2d2272b..7e4b63b65 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -3,6 +3,7 @@ module TrixiParticles using Reexport: @reexport using Adapt: Adapt +using Base: @propagate_inbounds using CSV: CSV using Dates using DataFrames: DataFrame diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 74a5f616d..2a3d45ea0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -23,15 +23,15 @@ difference of the coordinates, ``v_{ab} = v_a - v_b`` of the velocities of parti """ struct ContinuityDensity end -@inline function particle_density(v, system, particle) +@propagate_inbounds function particle_density(v, system, particle) particle_density(v, system.density_calculator, system, particle) end -@inline function particle_density(v, ::SummationDensity, system, particle) +@propagate_inbounds function particle_density(v, ::SummationDensity, system, particle) return system.cache.density[particle] end -@inline function particle_density(v, ::ContinuityDensity, system, particle) +@propagate_inbounds function particle_density(v, ::ContinuityDensity, system, particle) return v[end, particle] end diff --git a/src/general/system.jl b/src/general/system.jl index 674555aea..48d556eba 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -54,28 +54,31 @@ initialize!(system, neighborhood_search) = system @inline active_particles(system, ::Nothing) = eachparticle(system) # This should not be dispatched by system type. We always expect to get a column of `A`. -@inline function extract_svector(A, system, i) +@propagate_inbounds function extract_svector(A, system, i) extract_svector(A, Val(ndims(system)), i) end # Return the `i`-th column of the array `A` as an `SVector`. @inline function extract_svector(A, ::Val{NDIMS}, i) where {NDIMS} - return SVector(ntuple(@inline(dim->A[dim, i]), NDIMS)) + # Explicit bounds check, which can be removed by calling this function with `@inbounds` + @boundscheck checkbounds(A, NDIMS, i) + + # Assume inbounds access now + return SVector(ntuple(@inline(dim->@inbounds A[dim, i]), NDIMS)) end # Return `A[:, :, i]` as an `SMatrix`. @inline function extract_smatrix(A, system, particle) # Extract the matrix elements for this particle as a tuple to pass to SMatrix - return SMatrix{ndims(system), ndims(system)}( - # Convert linear index to Cartesian index - ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, - div(i - 1, ndims(system)) + 1, - particle]), - Val(ndims(system)^2))) + return SMatrix{ndims(system), + ndims(system)}(ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, + div(i - 1, ndims(system)) + 1, + particle]), + Val(ndims(system)^2))) end # Specifically get the current coordinates of a particle for all system types. -@inline function current_coords(u, system, particle) +@propagate_inbounds function current_coords(u, system, particle) return extract_svector(current_coordinates(u, system), system, particle) end @@ -91,7 +94,8 @@ end # This can be dispatched by system type. @inline initial_coordinates(system) = system.initial_condition.coordinates -@inline current_velocity(v, system, particle) = extract_svector(v, system, particle) +@propagate_inbounds current_velocity(v, system, particle) = extract_svector(v, system, + particle) @inline function current_acceleration(system, particle) # TODO: Return `dv` of solid particles diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index e0f07796c..fd0a9c3a5 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,7 +13,7 @@ function create_cache_density(ic, ::ContinuityDensity) return (;) end -@inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] +@propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] function write_u0!(u0, system::FluidSystem) (; initial_condition) = system diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 09bde77f9..95c9ff88e 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -1,9 +1,9 @@ # Unpack the neighboring systems viscosity to dispatch on the viscosity type -@inline function dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) viscosity = viscosity_model(particle_system, neighbor_system) return dv_viscosity(viscosity, particle_system, neighbor_system, @@ -12,10 +12,10 @@ sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) end -@inline function dv_viscosity(viscosity, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) return viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -105,12 +105,16 @@ function kinematic_viscosity(system, viscosity::ViscosityMorris) return viscosity.nu end -@inline function (viscosity::Union{ArtificialViscosityMonaghan, - ViscosityMorris})(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) +@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, + ViscosityMorris})(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, + pos_diff, distance, + sound_speed, + m_a, m_b, rho_a, rho_b, + grad_kernel) (; smoothing_length) = particle_system rho_mean = 0.5 * (rho_a + rho_b) @@ -250,4 +254,6 @@ function kinematic_viscosity(system, viscosity::ViscosityAdami) return viscosity.nu end -@inline viscous_velocity(v, system, particle) = current_velocity(v, system, particle) +@propagate_inbounds function viscous_velocity(v, system, particle) + return current_velocity(v, system, particle) +end diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c4ddc8b81..5d712235f 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -199,10 +199,10 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search return density_diffusion end -@inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, - v_particle_system, particle, neighbor, pos_diff, - distance, m_b, rho_a, rho_b, - particle_system::FluidSystem, grad_kernel) +@propagate_inbounds function density_diffusion!(dv, density_diffusion::DensityDiffusion, + v_particle_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system::FluidSystem, grad_kernel) # Density diffusion terms are all zero for distance zero distance < sqrt(eps()) && return diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 9e40eacfa..c74b25951 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -23,19 +23,23 @@ function interact!(dv, v_particle_system, u_particle_system, foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds particle_density(v_particle_system, particle_system, particle) + rho_b = @inbounds particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = 0.5 * (rho_a + rho_b) - # Determine correction values - viscosity_correction, pressure_correction, surface_tension_correction = free_surface_correction(correction, - particle_system, - rho_mean) + # Determine correction factors. + # This can be ignored, as these are all 1 when no correction is used. + (viscosity_correction, pressure_correction, + surface_tension_correction) = free_surface_correction(correction, particle_system, + rho_mean) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) # The following call is equivalent to # `p_a = particle_pressure(v_particle_system, particle_system, particle)` @@ -43,20 +47,23 @@ function interact!(dv, v_particle_system, u_particle_system, # Only when the neighbor system is a `BoundarySPHSystem` or a `TotalLagrangianSPHSystem` # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is # the pressure of the fluid particle. - p_a, p_b = particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) + p_a, p_b = @inbounds particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) dv_pressure = pressure_correction * pressure_acceleration(particle_system, neighbor_system, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) + # Propagate `@inbounds` to the viscosity function, which accesses particle data dv_viscosity_ = viscosity_correction * - dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + @inbounds dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, @@ -66,17 +73,19 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + - dv_adhesion[i] + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + + dv_surface_tension[i] + dv_adhesion[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end # TODO If variable smoothing_length is used, this should use the neighbor smoothing length - continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) + # Propagate `@inbounds` to the continuity equation, which accesses particle data + @inbounds continuity_equation!(dv, density_calculator, v_particle_system, + v_neighbor_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system, neighbor_system, grad_kernel) end # Debug example # periodic_box = neighborhood_search.periodic_box @@ -98,12 +107,12 @@ end end # This formulation was chosen to be consistent with the used pressure_acceleration formulations. -@inline function continuity_equation!(dv, density_calculator::ContinuityDensity, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::WeaklyCompressibleSPHSystem, - neighbor_system, grad_kernel) +@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::WeaklyCompressibleSPHSystem, + neighbor_system, grad_kernel) (; density_diffusion) = particle_system vdiff = current_velocity(v_particle_system, particle_system, particle) - @@ -121,9 +130,10 @@ end end end -@inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) +@propagate_inbounds function particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eb6c06c84..6d66ba06a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -211,7 +211,8 @@ end return ndims(system) + 1 end -@inline function particle_pressure(v, system::WeaklyCompressibleSPHSystem, particle) +@propagate_inbounds function particle_pressure(v, system::WeaklyCompressibleSPHSystem, + particle) return system.pressure[particle] end diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index e5f0c0e05..53cf860dd 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -131,7 +131,7 @@ Base.ntuple(f, ::Symbol) = ntuple(f, 2) # Make `extract_svector` work function TrixiParticles.current_coords(system::Val{:mock_system_tensor}, particle) - return TrixiParticles.extract_svector(current_coordinates[i], system, + return TrixiParticles.extract_svector(current_coordinates[i], Val(2), particle) end From 42e66b2a1afa25fbaa3f051216cfc7d9ec5b3679 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 Nov 2024 11:16:11 +0100 Subject: [PATCH 12/27] Clarify statement of need (#658) * update docs intro * improve text consistency * Update docs/src/index.md Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> * Update README.md Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> * Update README.md * Update README.md * Update README.md * Update index.md * Update README.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update docs/src/index.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update docs/src/index.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update docs/src/index.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update README.md * Update README.md --------- Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- README.md | 10 ++++++++-- docs/src/index.md | 10 +++++++++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 0fe4741b9..c90ed5b31 100644 --- a/README.md +++ b/README.md @@ -10,8 +10,14 @@ [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10797541.svg)](https://zenodo.org/doi/10.5281/zenodo.10797541) -**TrixiParticles.jl** is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in [Julia](https://julialang.org). -A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing. +**TrixiParticles.jl** is a high-performance numerical simulation framework for particle-based methods, focused on the simulation of complex multiphysics problems, and written in [Julia](https://julialang.org). + +TrixiParticles.jl focuses on the following use cases: +- Accurate and efficient physics-based modelling of complex multiphysics problems. +- Development of new particle-based methods and models. +- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work. + +It offers intuitive configuration, robust pre- and post-processing, and vendor-agnostic GPU-support based on the Julia package [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). [![YouTube](https://github.com/user-attachments/assets/dc2be627-a799-4bfd-9226-2077f737c4b0)](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s) diff --git a/docs/src/index.md b/docs/src/index.md index 42c2c42c6..27c1e9d47 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,14 @@ # TrixiParticles.jl -TrixiParticles.jl is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in Julia. A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing. Its features include: +**TrixiParticles.jl** is a high-performance particle simulation framework designed to overcome challenges of particle-based numerical methods in multiphysics applications. Existing frameworks often lack user-friendliness, involve complex configuration, and are not easily extensible for development of new methods. In the future we also want to provide seamless scalability from CPU to Exascale-level computing with GPU support. **TrixiParticles.jl** addresses these limitations with an intuitive interface, straightforward configuration, and an extensible design, facilitating efficient simulation setup and execution. + +TrixiParticles.jl focuses on the following use cases: + +- Development of new particle-based methods and models. By providing an extensible architecture to incorporate additional particle methods easily and not focusing on a single model or numerical method. +- Accurate, reliable and efficient physics-based modelling of complex multiphysics problems by providing a flexible configuration system, tools, high performance and a wide range of validation and test cases. +- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work through extensive documentation, community engagement and readable configuration files. + +Its features include: ## Features - Incompressible Navier-Stokes From d52b1d51a46581fdca4a94e2fa368b65fb875fea Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 18 Nov 2024 14:41:30 +0100 Subject: [PATCH 13/27] Make simulations run with `Float32` (#662) * Replace hardcoded Float64 values * Avoid Float64 `eps()` * Remove the last Float64 occurrence in fluid-fluid kernel * Fix tests * Make all smoothing kernels preserve types and add test for this * Make eps(typeof(h)) consistent * Reformat code * Fix tests --- src/general/corrections.jl | 4 +- src/general/smoothing_kernels.jl | 75 ++++++++++++------- src/schemes/fluid/viscosity.jl | 2 +- .../density_diffusion.jl | 4 +- .../fluid/weakly_compressible_sph/rhs.jl | 2 +- test/general/smoothing_kernels.jl | 35 ++++++++- .../schemes/solid/total_lagrangian_sph/rhs.jl | 1 + test/systems/solid_system.jl | 3 +- 8 files changed, 89 insertions(+), 37 deletions(-) diff --git a/src/general/corrections.jl b/src/general/corrections.jl index cf6c8bed6..b1c66a06a 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -40,11 +40,11 @@ end k = correction.rho0 / rho_mean # Viscosity, pressure, surface_tension - return k, 1.0, k + return k, 1, k end @inline function free_surface_correction(correction, particle_system, rho_mean) - return 1.0, 1.0, 1.0 + return 1, 1, 1 end @doc raw""" diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index 20268ebad..862e58ba1 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -3,7 +3,8 @@ abstract type SmoothingKernel{NDIMS} end @inline Base.ndims(::SmoothingKernel{NDIMS}) where {NDIMS} = NDIMS @inline function kernel_grad(kernel, pos_diff, distance, h) - distance < sqrt(eps()) && return zero(pos_diff) + # TODO Use `eps` relative to `h` to allow scaling of simulations + distance < sqrt(eps(typeof(h))) && return zero(pos_diff) return kernel_deriv(kernel, distance, h) / distance * pos_diff end @@ -98,7 +99,8 @@ end @inline compact_support(::GaussianKernel, h) = 3 * h @inline normalization_factor(::GaussianKernel{2}, h) = 1 / (pi * h^2) -@inline normalization_factor(::GaussianKernel{3}, h) = 1 / (pi^(3 / 2) * h^3) +# First convert `pi` to the type of `h` to preserve the type of `h` +@inline normalization_factor(::GaussianKernel{3}, h) = 1 / (oftype(h, pi)^(3 // 2) * h^3) @doc raw""" SchoenbergCubicSplineKernel{NDIMS}() @@ -136,8 +138,9 @@ struct SchoenbergCubicSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end @muladd @inline function kernel(kernel::SchoenbergCubicSplineKernel, r::Real, h) q = r / h - # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl - result = 1 / 4 * (2 - q)^3 + # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl. + # Use `//` to preserve the type of `q`. + result = 1 // 4 * (2 - q)^3 result = result - (q < 1) * (1 - q)^3 # Zero out result if q >= 2 @@ -151,7 +154,8 @@ end q = r * inner_deriv # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl - result = -3 / 4 * (2 - q)^2 + # Use `//` to preserve the type of `q`. + result = -3 // 4 * (2 - q)^2 result = result + 3 * (q < 1) * (1 - q)^2 # Zero out result if q >= 2 @@ -164,7 +168,8 @@ end @inline compact_support(::SchoenbergCubicSplineKernel, h) = 2 * h @inline normalization_factor(::SchoenbergCubicSplineKernel{1}, h) = 2 / 3h -@inline normalization_factor(::SchoenbergCubicSplineKernel{2}, h) = 10 / (7 * pi * h^2) +# `7 * pi` is always `Float64`. `pi * h^2 * 7` preserves the type of `h`. +@inline normalization_factor(::SchoenbergCubicSplineKernel{2}, h) = 10 / (pi * h^2 * 7) @inline normalization_factor(::SchoenbergCubicSplineKernel{3}, h) = 1 / (pi * h^3) @doc raw""" @@ -213,17 +218,19 @@ struct SchoenbergQuarticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end # is a higher priority than exact precision. @fastpow @muladd @inline function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h) q = r / h - q5_2 = (5 / 2 - q) - q3_2 = (3 / 2 - q) - q1_2 = (1 / 2 - q) + + # Preserve the type of `q` + q5_2 = (5 // 2 - q) + q3_2 = (3 // 2 - q) + q1_2 = (1 // 2 - q) # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl result = q5_2^4 - result = result - 5 * (q < 3 / 2) * q3_2^4 - result = result + 10 * (q < 1 / 2) * q1_2^4 + result = result - 5 * (q < 3 // 2) * q3_2^4 + result = result + 10 * (q < 1 // 2) * q1_2^4 # Zero out result if q >= 5/2 - result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result, zero(result)) + result = ifelse(q < 5 // 2, normalization_factor(kernel, h) * result, zero(result)) return result end @@ -232,17 +239,19 @@ end r::Real, h) inner_deriv = 1 / h q = r * inner_deriv - q5_2 = 5 / 2 - q - q3_2 = 3 / 2 - q - q1_2 = 1 / 2 - q + + # Preserve the type of `q` + q5_2 = 5 // 2 - q + q3_2 = 3 // 2 - q + q1_2 = 1 // 2 - q # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl result = -4 * q5_2^3 - result = result + 20 * (q < 3 / 2) * q3_2^3 - result = result - 40 * (q < 1 / 2) * q1_2^3 + result = result + 20 * (q < 3 // 2) * q3_2^3 + result = result - 40 * (q < 1 // 2) * q1_2^3 # Zero out result if q >= 5/2 - result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result * inner_deriv, + result = ifelse(q < 5 // 2, normalization_factor(kernel, h) * result * inner_deriv, zero(result)) return result @@ -251,8 +260,9 @@ end @inline compact_support(::SchoenbergQuarticSplineKernel, h) = 2.5 * h @inline normalization_factor(::SchoenbergQuarticSplineKernel{1}, h) = 1 / 24h -@inline normalization_factor(::SchoenbergQuarticSplineKernel{2}, h) = 96 / (1199 * pi * h^2) -@inline normalization_factor(::SchoenbergQuarticSplineKernel{3}, h) = 1 / (20 * pi * h^3) +# `1199 * pi` is always `Float64`. `pi * h^2 * 1199` preserves the type of `h`. +@inline normalization_factor(::SchoenbergQuarticSplineKernel{2}, h) = 96 / (pi * h^2 * 1199) +@inline normalization_factor(::SchoenbergQuarticSplineKernel{3}, h) = 1 / (pi * h^3 * 20) @doc raw""" SchoenbergQuinticSplineKernel{NDIMS}() @@ -329,8 +339,9 @@ end @inline compact_support(::SchoenbergQuinticSplineKernel, h) = 3 * h @inline normalization_factor(::SchoenbergQuinticSplineKernel{1}, h) = 1 / 120h -@inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (478 * pi * h^2) -@inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (120 * pi * h^3) +# `478 * pi` is always `Float64`. `pi * h^2 * 478` preserves the type of `h`. +@inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (pi * h^2 * 478) +@inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (pi * h^3 * 120) abstract type WendlandKernel{NDIMS} <: SmoothingKernel{NDIMS} end @@ -400,7 +411,8 @@ end end @inline normalization_factor(::WendlandC2Kernel{2}, h) = 7 / (pi * h^2) -@inline normalization_factor(::WendlandC2Kernel{3}, h) = 21 / (2pi * h^3) +# `2 * pi` is always `Float64`. `pi * h^3 * 2` preserves the type of `h`. +@inline normalization_factor(::WendlandC2Kernel{3}, h) = 21 / (pi * h^3 * 2) @doc raw""" WendlandC4Kernel{NDIMS}() @@ -449,8 +461,10 @@ end @fastpow @muladd @inline function kernel_deriv(kernel::WendlandC4Kernel, r::Real, h) q = r / h - term1 = (1 - q)^6 * (6 + 70 / 3 * q) - term2 = 6 * (1 - q)^5 * (1 + 6q + 35 / 3 * q^2) + + # Use `//` to preserve the type of `q` + term1 = (1 - q)^6 * (6 + 70 // 3 * q) + term2 = 6 * (1 - q)^5 * (1 + 6q + 35 // 3 * q^2) derivative = term1 - term2 # Zero out result if q >= 1 @@ -461,7 +475,8 @@ end end @inline normalization_factor(::WendlandC4Kernel{2}, h) = 9 / (pi * h^2) -@inline normalization_factor(::WendlandC4Kernel{3}, h) = 495 / (32pi * h^3) +# `32 * pi` is always `Float64`. `pi * h^2 * 32` preserves the type of `h`. +@inline normalization_factor(::WendlandC4Kernel{3}, h) = 495 / (pi * h^3 * 32) @doc raw""" WendlandC6Kernel{NDIMS}() @@ -521,8 +536,9 @@ end return result end -@inline normalization_factor(::WendlandC6Kernel{2}, h) = 78 / (7pi * h^2) -@inline normalization_factor(::WendlandC6Kernel{3}, h) = 1365 / (64pi * h^3) +# `7 * pi` is always `Float64`. `pi * h^2 * 7` preserves the type of `h`. +@inline normalization_factor(::WendlandC6Kernel{2}, h) = 78 / (pi * h^2 * 7) +@inline normalization_factor(::WendlandC6Kernel{3}, h) = 1365 / (pi * h^3 * 64) @doc raw""" Poly6Kernel{NDIMS}() @@ -588,7 +604,8 @@ end @inline compact_support(::Poly6Kernel, h) = h @inline normalization_factor(::Poly6Kernel{2}, h) = 4 / (pi * h^2) -@inline normalization_factor(::Poly6Kernel{3}, h) = 315 / (64pi * h^3) +# `64 * pi` is always `Float64`. `pi * h^3 * 64` preserves the type of `h`. +@inline normalization_factor(::Poly6Kernel{3}, h) = 315 / (pi * h^3 * 64) @doc raw""" SpikyKernel{NDIMS}() diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 95c9ff88e..1a0665a71 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -117,7 +117,7 @@ end grad_kernel) (; smoothing_length) = particle_system - rho_mean = 0.5 * (rho_a + rho_b) + rho_mean = (rho_a + rho_b) / 2 v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index 5d712235f..c10740ab8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -178,7 +178,7 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search foreach_point_neighbor(system, system, system_coords, system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return + distance < sqrt(eps(typeof(distance))) && return rho_a = particle_density(v, system, particle) rho_b = particle_density(v, system, neighbor) @@ -204,7 +204,7 @@ end pos_diff, distance, m_b, rho_a, rho_b, particle_system::FluidSystem, grad_kernel) # Density diffusion terms are all zero for distance zero - distance < sqrt(eps()) && return + distance < sqrt(eps(typeof(distance))) && return (; delta) = density_diffusion (; smoothing_length, state_equation) = particle_system diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index c74b25951..0fe58ab21 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -28,7 +28,7 @@ function interact!(dv, v_particle_system, u_particle_system, # in this hot loop to avoid bounds checking when extracting particle quantities. rho_a = @inbounds particle_density(v_particle_system, particle_system, particle) rho_b = @inbounds particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = 0.5 * (rho_a + rho_b) + rho_mean = (rho_a + rho_b) / 2 # Determine correction factors. # This can be ignored, as these are all 1 when no correction is used. diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index 5d6f15c26..cb9c98f4b 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -90,4 +90,37 @@ end end end -end + + @testset verbose=false "Return Type" begin + # Test that the return type of the kernel and kernel derivative preserve + # the input type. We don't want to return `Float64` when working with `Float32`. + kernels = [ + GaussianKernel, + SchoenbergCubicSplineKernel, + SchoenbergQuarticSplineKernel, + SchoenbergQuinticSplineKernel, + WendlandC2Kernel, + WendlandC4Kernel, + WendlandC6Kernel, + SpikyKernel, + Poly6Kernel + ] + + # Test different smoothing length types + smoothing_lengths = (0.5, 0.5f0) + + @testset "$kernel_type" for kernel_type in kernels + for ndims in 2:3 + kernel_ = kernel_type{ndims}() + + for h in smoothing_lengths + result = TrixiParticles.kernel(kernel_, h / 2, h) + @test typeof(result) == typeof(h) + + result = TrixiParticles.kernel_deriv(kernel_, h / 2, h) + @test typeof(result) == typeof(h) + end + end + end + end +end; diff --git a/test/schemes/solid/total_lagrangian_sph/rhs.jl b/test/schemes/solid/total_lagrangian_sph/rhs.jl index ea2474d89..b69f835ca 100644 --- a/test/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/test/schemes/solid/total_lagrangian_sph/rhs.jl @@ -114,6 +114,7 @@ nothing end TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv + Base.eps(::Type{Val{:mock_smoothing_length}}) = eps() #### Verification systems = [system, system_gpu] diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index 53cf860dd..18b887b49 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -171,6 +171,7 @@ function TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) return kernel_derivative end + Base.eps(::Type{Val{:mock_smoothing_length}}) = eps() # Compute deformation gradient deformation_grad = ones(2, 2, 2) @@ -371,4 +372,4 @@ @test isapprox(von_mises_stress[1], 1.4257267477533202, atol=1e-14) @test isapprox(reference_stress_tensor, cauchy_stress, atol=1e-6) end -end +end; From b8eb490788977d6e294f6c3f811228e437c18a25 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 Nov 2024 16:38:29 +0100 Subject: [PATCH 14/27] Fixes for Open Boundary Systems for free surfaces (#609) * set test up for 1.11 * fix * typo * fix again * remove soundspeed from OBS * skip empty system * fix test * fix tests * fix tests * fix bug * fix * check dimensionality of reference functions * propagate characteristics * update * cleanup * update * Increase errors for 1.11 * Fix invalidations * Fix tests * Update ci.yml * revert * Update ci.yml * Update test/validation/validation.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * revert changes that had no benefit * update * cleanup * include in test run * remove redundancy * revert * fix tests * fix * fix test * fix test * fix the test * fix error * Update src/schemes/boundary/open_boundary/method_of_characteristics.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update test/schemes/boundary/open_boundary/boundary_zone.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update src/setups/extrude_geometry.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update src/schemes/boundary/open_boundary/method_of_characteristics.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update src/schemes/boundary/open_boundary/boundary_zones.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update examples/fluid/pipe_flow_3d.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * fix test * fix test * format * fix test * format --------- Co-authored-by: LasNikas Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- examples/fluid/pipe_flow_2d.jl | 15 ++++--- examples/fluid/pipe_flow_3d.jl | 45 +++++++++++++++++++ .../boundary/open_boundary/boundary_zones.jl | 17 +++---- .../method_of_characteristics.jl | 23 +++++++--- src/schemes/boundary/open_boundary/system.jl | 28 ++++++++++-- src/setups/extrude_geometry.jl | 4 +- test/examples/examples.jl | 8 ++++ .../boundary/open_boundary/boundary_zone.jl | 4 +- 8 files changed, 116 insertions(+), 28 deletions(-) create mode 100644 examples/fluid/pipe_flow_3d.jl diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index af558133f..c8778f485 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -48,7 +48,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi # Shift pipe walls in negative x-direction for the inflow pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers -n_buffer_particles = 4 * pipe.n_particles_per_dimension[2] +n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1) # ========================================================================================== # ==== Fluid @@ -77,14 +77,16 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi # ========================================================================================== # ==== Open Boundary -function velocity_function(pos, t) + +function velocity_function2d(pos, t) # Use this for a time-dependent inflow velocity # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) return SVector(prescribed_velocity, 0.0) end -inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, +plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = InFlow(; plane=plane_in, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, @@ -92,9 +94,10 @@ open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) -outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), +plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) +outflow = OutFlow(; plane=plane_out, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) @@ -103,7 +106,7 @@ open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) # ========================================================================================== # ==== Boundary diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl new file mode 100644 index 000000000..fb2cc7f0a --- /dev/null +++ b/examples/fluid/pipe_flow_3d.jl @@ -0,0 +1,45 @@ +# 3D channel flow simulation with open boundaries. + +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +function velocity_function3d(pos, t) + # Use this for a time-dependent inflow velocity + # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) + + return SVector(prescribed_velocity, 0.0, 0.0) +end + +domain_size = (1.0, 0.4, 0.4) + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) + +flow_direction = [1.0, 0.0, 0.0] + +# setup simulation +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + domain_size=domain_size, boundary_size=boundary_size, + flow_direction=flow_direction, faces=(false, false, true, true, true, true), + tspan=tspan, smoothing_kernel=WendlandC2Kernel{3}(), + reference_velocity=velocity_function3d, + plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0], + [0.0, 0.0, domain_size[3]]), + plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0], + [domain_size[1], 0.0, domain_size[3]])) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index cc9a77157..ad72045d6 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -116,11 +116,11 @@ struct InFlow{NDIMS, IC, S, ZO, ZW, FD} if !isapprox(abs(dot_), 1.0, atol=1e-7) throw(ArgumentError("`flow_direction` is not normal to inflow plane")) - else - # Flip the normal vector to point in the opposite direction of `flow_direction` - spanning_set[:, 1] .*= -sign(dot_) end + # Flip the normal vector to point in the opposite direction of `flow_direction` + spanning_set[:, 1] .*= -sign(dot_) + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) # Remove particles outside the boundary zone. @@ -251,11 +251,11 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} if !isapprox(abs(dot_), 1.0, atol=1e-7) throw(ArgumentError("`flow_direction` is not normal to outflow plane")) - else - # Flip the normal vector to point in `flow_direction` - spanning_set[:, 1] .*= sign(dot_) end + # Flip the normal vector to point in `flow_direction` + spanning_set[:, 1] .*= sign(dot_) + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) # Remove particles outside the boundary zone. @@ -289,8 +289,9 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) edge1 = plane_points[2] - plane_points[1] edge2 = plane_points[3] - plane_points[1] - if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) - throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal")) + # Check if the edges are linearly dependent (to avoid degenerate planes) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("the vectors `AB` and `AC` must not be collinear")) end # Calculate normal vector of plane diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 38fdd477a..f803bf9ca 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -8,8 +8,9 @@ about the method see [description below](@ref method_of_characteristics). """ struct BoundaryModelLastiwka end -@inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka, - v, u, v_ode, u_ode, semi, t) +# Called from update callback via `update_open_boundary_eachstep!` +@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, + v, u, v_ode, u_ode, semi, t) (; density, pressure, cache, flow_direction, reference_velocity, reference_pressure, reference_density) = system @@ -45,10 +46,13 @@ struct BoundaryModelLastiwka end return system end +# Called from semidiscretization function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t) @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end + + return system end # ==== Characteristics @@ -98,9 +102,16 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end end - characteristics[1, particle] = avg_J1 / counter - characteristics[2, particle] = avg_J2 / counter - characteristics[3, particle] = avg_J3 / counter + # To prevent NANs here if the boundary has not been in contact before. + if counter > 0 + characteristics[1, particle] = avg_J1 / counter + characteristics[2, particle] = avg_J2 / counter + characteristics[3, particle] = avg_J3 / counter + else + characteristics[1, particle] = 0 + characteristics[2, particle] = 0 + characteristics[3, particle] = 0 + end else characteristics[1, particle] /= volume[particle] characteristics[2, particle] /= volume[particle] @@ -164,7 +175,7 @@ end @inline function prescribe_conditions!(characteristics, particle, ::OutFlow) # J3 is prescribed (i.e. determined from the exterior of the domain). - # J1 and J2 is transimtted from the domain interior. + # J1 and J2 is transmitted from the domain interior. characteristics[3, particle] = zero(eltype(characteristics)) return characteristics diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9d2c79b00..c8d80ca8f 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -77,6 +77,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "an array where the ``i``-th column holds the velocity of particle ``i`` " * "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) else + if reference_velocity isa Function + test_result = reference_velocity(zeros(NDIMS), 0.0) + if length(test_result) != NDIMS + throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) + end + end reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) end @@ -86,6 +92,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its pressure, " * "a vector holding the pressure of each particle, or a scalar")) else + if reference_pressure isa Function + test_result = reference_pressure(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_pressure` function must be a scalar function")) + end + end reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) end @@ -95,6 +107,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its density, " * "a vector holding the density of each particle, or a scalar")) else + if reference_density isa Function + test_result = reference_density(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_density` function must be a scalar function")) + end + end reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end @@ -190,10 +208,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ # Update density, pressure and velocity based on the characteristic variables. # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 - @trixi_timeit timer() "update quantities" update_quantities!(system, - system.boundary_model, - v, u, v_ode, u_ode, - semi, t) + @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, + system.boundary_model, + v, u, + v_ode, + u_ode, + semi, t) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index c666ec3d7..9b31c463c 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -177,8 +177,8 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing) edge2 = point3_ - point1_ # Check if the points are collinear - if norm(cross(edge1, edge2)) == 0 - throw(ArgumentError("the points must not be collinear")) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("the vectors `AB` and `AC` must not be collinear")) end # Determine the number of points along each edge diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 2e5462a35..3673e5bd5 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -238,6 +238,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/pipe_flow_3d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "pipe_flow_3d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/lid_driven_cavity_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index ad1b69c4b..2ed6e2810 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -187,10 +187,10 @@ end @testset verbose=true "Illegal Inputs" begin - no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.4, 0.9, -0.15]] flow_direction = [0.0, 0.0, 1.0] - error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal" + error_str = "the vectors `AB` and `AC` must not be collinear" @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, particle_spacing=0.1, From 7df5b1821d365b4ae64e9339d86965e40d0a01ba Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 17:34:41 +0100 Subject: [PATCH 15/27] add callback --- src/callbacks/callbacks.jl | 1 + src/callbacks/refinement.jl | 123 ++++++++++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 src/callbacks/refinement.jl 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..03f7f6fe7 --- /dev/null +++ b/src/callbacks/refinement.jl @@ -0,0 +1,123 @@ +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_update!, + save_positions=(false, false)) + else + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(refinement_callback, refinement_callback, + initialize=initial_update!, + 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) + 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 + + # Update NHS + @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi) + + # 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 From 022c6df6c31cbfae09ab143a5ffc4b4251299cdc Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 12:38:05 +0100 Subject: [PATCH 16/27] First prototype of variable smoothing length --- src/general/corrections.jl | 6 +-- src/general/density_calculators.jl | 2 +- src/general/system.jl | 25 +++++------- src/schemes/boundary/rhs.jl | 2 +- .../fluid/entropically_damped_sph/rhs.jl | 9 ++--- src/schemes/fluid/surface_tension.jl | 6 +-- src/schemes/fluid/viscosity.jl | 38 ++++++++++++------- .../density_diffusion.jl | 13 ++++--- .../fluid/weakly_compressible_sph/system.jl | 31 +++++++++++++-- src/schemes/solid/total_lagrangian_sph/rhs.jl | 4 +- .../solid/total_lagrangian_sph/system.jl | 2 +- 11 files changed, 84 insertions(+), 54 deletions(-) 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/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/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/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3644c3630..862471ed1 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,7 +30,7 @@ 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, m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, @@ -72,8 +72,6 @@ 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 - volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -88,8 +86,9 @@ end 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 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/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1a0665a71..aa8bc2c0f 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,18 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system - 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 +242,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 +262,7 @@ end return visc .* v_diff end -function kinematic_viscosity(system, viscosity::ViscosityAdami) +function kinematic_viscosity(system, viscosity::ViscosityAdami, particle) 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/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 6d66ba06a..df1d65d5d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -45,14 +45,13 @@ 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} + 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,6 +117,9 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, n_particles)..., + cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, @@ -123,7 +127,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, 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) = (;) @@ -155,6 +160,26 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end +function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length) +end + +function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length=smoothing_length * ones(n_particles)) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..ed1f4f745 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. 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' From fe4934aaf77dcbd266327765189722d456dc92c3 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 17/27] Make it work with EDAC --- src/general/semidiscretization.jl | 8 +++++--- .../boundary/dummy_particles/dummy_particles.jl | 6 +++++- src/schemes/fluid/entropically_damped_sph/system.jl | 13 ++++++++++++- src/schemes/fluid/weakly_compressible_sph/system.jl | 2 +- src/visualization/write2vtk.jl | 2 +- 5 files changed, 24 insertions(+), 7 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index cf9df1b5e..f734a6191 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -137,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) @@ -410,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/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 8b7e507ed..345e3c807 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -69,6 +69,10 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, smoothing_length, viscosity, correction, cache) end +function smoothing_length(boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + @doc raw""" AdamiPressureExtrapolation(; pressure_offset=0.0) @@ -480,7 +484,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 + diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 94064555d..304baaaa2 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -49,7 +49,6 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K - smoothing_length :: ELTYPE sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -146,6 +145,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index df1d65d5d..5a4a0d103 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -123,7 +123,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, 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 From 17f070d981648fcc316691463cd28fd15432cdfb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 18/27] add h-factor --- .../fluid/entropically_damped_sph/rhs.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 23 +++++---------- src/schemes/fluid/fluid.jl | 22 +++++++++++++++ .../fluid/weakly_compressible_sph/system.jl | 28 +++---------------- 4 files changed, 34 insertions(+), 41 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 862471ed1..71a639974 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -87,7 +87,7 @@ end eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) smoothing_length_average = 0.5 * (smoothing_length(particle_system, particle) + - 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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 304baaaa2..485078c9e 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -101,12 +101,15 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, 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...) return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, - smoothing_kernel, smoothing_length, sound_speed, - viscosity, nu_edac, acceleration_, nothing, - pressure_acceleration, transport_velocity, - source_terms, buffer, particle_refinement, cache) + 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) @@ -145,18 +148,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fd0a9c3a5..1b0413010 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,8 +13,30 @@ 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 + return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)), + smoothing_length_factor=smoothng_length_factor) +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 diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 5a4a0d103..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,8 +44,8 @@ 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, PR, 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} @@ -118,8 +118,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) cache = (; - create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, n_particles)..., - cache...) + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, @@ -160,26 +160,6 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end -function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length) -end - -function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length=smoothing_length * ones(n_particles)) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time From b6a51081c527e2b808b3143f16f795e29b68e3b7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 13:04:47 +0100 Subject: [PATCH 19/27] adapt rhs --- .../fluid/entropically_damped_sph/rhs.jl | 210 ++++++++++++++++-- src/schemes/fluid/fluid.jl | 7 +- src/schemes/fluid/pressure_acceleration.jl | 4 +- src/schemes/fluid/transport_velocity.jl | 25 ++- .../fluid/weakly_compressible_sph/rhs.jl | 3 +- src/schemes/solid/total_lagrangian_sph/rhs.jl | 3 +- 6 files changed, 230 insertions(+), 22 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 71a639974..481030ec9 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -32,7 +32,8 @@ function interact!(dv, v_particle_system, u_particle_system, 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,28 @@ 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_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,9 +78,45 @@ 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) +@inline function pressure_evolution!(dv, particle_system, 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_system.particle_refinement, + particle) + beta_inv_b = beta_correction(neighbor_system, neighbor_system.particle_refinement, + 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) + + 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 + +function beta_correction(particle_system, ::Nothing, particle) + return one(eltype(particle_system)) +end + +function beta_correction(particle_system, ::ParticleRefinement, 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_b, rho_b) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -79,9 +124,6 @@ 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) @@ -101,11 +143,109 @@ 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, ::ParticleRefinement, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_b, rho_b) + # EDAC pressure evolution + pressure_diff = p_a - p_b - return dv + # 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 * particle_system.smoothing_length[particle] / tmp + nu_edac_a = alpha_edac * sound_speed * particle_system.smoothing_length[neighbor] / tmp + + nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) + + # TODO: Use wrapped version + grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[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_damping_term(particle_system, neighbor_system, ::ParticleRefinement, + particle, neighbor, pos_diff, distance, 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, ::ParticleRefinement, + 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_refinement, particle) + beta_inv_b = beta_correction(particle_system, particle_refinement, 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 = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[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) + 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(eltype(particle_system)) +end + +@inline function velocity_correction(particle_system, neighbor_system, + particle_refinement::ParticleRefinement, + 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) + + 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_refinement, 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,:]`. @@ -124,3 +264,45 @@ 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::ParticleRefinement, + 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_refinement, particle) * + (p_a - p_a_avg) / rho_a^2 + P_b = beta_correction(particle_system, particle_refinement, neighbor) * + (p_b - p_b_avg) / rho_b^2 + + # TODO: Use wrapped version + grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[neighbor]) + + return -m_b * (P_a * grad_kernel_a + P_b * grad_kernel_b) +end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 1b0413010..9b8896c3e 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -111,7 +111,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) @@ -120,9 +120,10 @@ 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) + pos_diff, distance, v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) 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/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bc8fdf005..90b744f1d 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, + ::ParticleRefinement, 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(particle_system, particle_refinement, particle) / rho_a + factor_b = beta_correction(neighbor_system, particle_refinement, 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, 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/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index ed1f4f745..114f72aed 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -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) From 1e2651893aa4301951ec944fa8744fc395575383 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 13:07:20 +0100 Subject: [PATCH 20/27] =?UTF-8?q?ada=C3=BCt=20viscosity?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../dummy_particles/dummy_particles.jl | 11 +-- .../fluid/entropically_damped_sph/rhs.jl | 81 ++++++++----------- src/schemes/fluid/fluid.jl | 5 +- src/schemes/fluid/transport_velocity.jl | 6 +- src/schemes/fluid/viscosity.jl | 41 ++++++++-- 5 files changed, 81 insertions(+), 63 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 345e3c807..ece912eeb 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -585,12 +585,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/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 481030ec9..b075813f2 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -64,7 +64,8 @@ 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) - pressure_evolution!(dv, particle_system, neighbor_system, v_diff, grad_kernel, + 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) @@ -79,6 +80,7 @@ function interact!(dv, v_particle_system, u_particle_system, end @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) @@ -87,16 +89,14 @@ end # 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_system.particle_refinement, - particle) - beta_inv_b = beta_correction(neighbor_system, neighbor_system.particle_refinement, - neighbor) + 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) + + 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, @@ -105,18 +105,21 @@ end return dv end +@inline beta_correction(system, particle) = one(eltype(system)) -function beta_correction(particle_system, ::Nothing, particle) - return one(eltype(particle_system)) +@inline function beta_correction(system::FluidSystem, particle) + beta_correction(system, system.particle_refinement, particle) end -function beta_correction(particle_system, ::ParticleRefinement, particle) +@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_b, rho_b) + 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 @@ -146,9 +149,9 @@ function pressure_damping_term(particle_system, neighbor_system, ::Nothing, return volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) end -function pressure_damping_term(particle_system, neighbor_system, ::ParticleRefinement, +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_b, rho_b) + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel_a) # EDAC pressure evolution pressure_diff = p_a - p_b @@ -158,24 +161,17 @@ function pressure_damping_term(particle_system, neighbor_system, ::ParticleRefin # 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 * particle_system.smoothing_length[particle] / tmp - nu_edac_a = alpha_edac * sound_speed * particle_system.smoothing_length[neighbor] / tmp + nu_edac_a = alpha_edac * sound_speed * smoothing_length(particle_system, particle) / tmp + nu_edac_a = alpha_edac * sound_speed * smoothing_length(neighbor_system, particle) / tmp nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) - # TODO: Use wrapped version - grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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_damping_term(particle_system, neighbor_system, ::ParticleRefinement, - particle, neighbor, pos_diff, distance, m_b, rho_b) -end function pressure_reduction(particle_system, neighbor_system, ::Nothing, v_particle_system, v_neighbor_system, @@ -184,12 +180,12 @@ function pressure_reduction(particle_system, neighbor_system, ::Nothing, return zero(eltype(particle_system)) end -function pressure_reduction(particle_system, neighbor_system, ::ParticleRefinement, +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_refinement, particle) - beta_inv_b = beta_correction(particle_system, particle_refinement, neighbor) + 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) @@ -197,10 +193,8 @@ function pressure_reduction(particle_system, neighbor_system, ::ParticleRefineme 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 = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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) @@ -224,12 +218,11 @@ end v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) - return zero(eltype(particle_system)) + return zero(pos_diff) end @inline function velocity_correction(particle_system, neighbor_system, - particle_refinement::ParticleRefinement, - pos_diff, distance, + particle_refinement, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) @@ -242,7 +235,7 @@ end 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_refinement, particle) + 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 @@ -286,23 +279,17 @@ end end @inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, - particle_refinement::ParticleRefinement, - neighbor_system, particle, neighbor, m_a, m_b, p_a, - p_b, - rho_a, rho_b, pos_diff, distance, W_a, correction) + 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_refinement, particle) * - (p_a - p_a_avg) / rho_a^2 - P_b = beta_correction(particle_system, particle_refinement, neighbor) * - (p_b - p_b_avg) / rho_b^2 + 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 - # TODO: Use wrapped version - grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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/fluid.jl b/src/schemes/fluid/fluid.jl index 9b8896c3e..f287cabba 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -124,6 +124,7 @@ end 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, - pos_diff, distance, 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 diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 90b744f1d..5270d1210 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -103,7 +103,7 @@ end end @inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, - ::ParticleRefinement, pos_diff, distance, + 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) @@ -112,8 +112,8 @@ end 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(particle_system, particle_refinement, particle) / rho_a - factor_b = beta_correction(neighbor_system, particle_refinement, neighbor) / rho_b + factor_a = beta_correction(particle_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)' diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index aa8bc2c0f..222fa6fd8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -112,12 +112,9 @@ end v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, - m_a, m_b, rho_a, rho_b, - grad_kernel) - (; smoothing_length) = particle_system - - rho_mean = (rho_a + rho_b) / 2 + sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + rho_mean = 0.5 * (rho_a + rho_b) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -220,6 +217,38 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + 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)) + + 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) From a5582f94b270049c35615c67217ac882d71cacbb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:06:13 +0100 Subject: [PATCH 21/27] add `ParticleRefinement` --- src/TrixiParticles.jl | 3 +- src/refinement/refinement.jl | 136 ++++++++++++++++++++++++-- src/refinement/refinement_criteria.jl | 64 ++++++++++++ 3 files changed, 196 insertions(+), 7 deletions(-) create mode 100644 src/refinement/refinement_criteria.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 7e4b63b65..d87a03796 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -45,6 +45,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. @@ -73,7 +74,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/refinement/refinement.jl b/src/refinement/refinement.jl index bf886ac12..ae435e481 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,9 +1,133 @@ -struct ParticleRefinement - n_particles_before_resize :: Ref{Int} - n_new_particles :: Ref{Int} - delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles +include("refinement_criteria.jl") + +struct ParticleRefinement{RC, ELTYPE} + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles +end + +function ParticleRefinement(; max_spacing_ratio, + refinement_criteria=SpatialRefinementCriterion()) + mass_ref = Vector{eltype(max_spacing_ratio)}() + + if refinement_criteria isa Tuple + refinement_criteria = (refinement_criteria,) + end + return ParticleRefinement(refinement_criteria, max_spacing_ratio, mass_ref) +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)) + + return system +end + +function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + check_refinement_criteria!(semi, v_ode, u_ode) + + # Update the spacing of particles (Algorthm 1) + update_particle_spacing(semi, u_ode) + + # Split the particles (Algorithm 2) + + # Merge the particles (Algorithm 3) + + # Shift the particles + + # Correct the particles + + # Update smoothing lengths + + # Resize neighborhood search + + return semi end -function ParticleRefinement() - return ParticleRefinement(Ref(0), Ref(0), Bool[]) +function update_particle_spacing(semi, u_ode) + foreach_system(semi) do system + u = wrap_u(u_ode, system, semi) + update_particle_spacing(system, u, semi) + end +end + +# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl` +@inline update_particle_spacing(system, u, semi) = systemerror + +@inline function update_particle_spacing(system::FluidSystem, u, semi) + update_particle_spacing(system, system.particle_refinement, u, semi) +end + +@inline update_particle_spacing(system, ::Nothing, u, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, u, semi) + (; smoothing_length, smoothing_length_factor) = system.cache + (; mass_ref) = particle_refinement + + system_coords = current_coordinates(u, system) + + for particle in eachparticle(system) + dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, 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] = system.density[particle] * new_spacing^(ndims(system)) + end + + return system +end + +@inline function min_max_avg_spacing(system, semi, system_coords, particle) + dp_min = Inf + dp_max = zero(eltype(system)) + dp_avg = zero(eltype(system)) + counter_neighbors = 0 + + foreach_system(semi) do neighbor_system + neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, + semi) + + u_neighbor_system = wrap_u(u_ode, neighbor, 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 + + 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..de4bab43a --- /dev/null +++ b/src/refinement/refinement_criteria.jl @@ -0,0 +1,64 @@ +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, 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_neighbor = particle_spacing(neighbor_system, neighbor) + dp_particle = smoothing_length[particle] / smoothing_length_factor + + smoothing_length[particle] = smoothing_length_factor * min(dp_neighbor, dp_particle) + end + + return particle_system +end From 6f97bf961265c328e6d38c16a99325c73e3506db Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:08:11 +0100 Subject: [PATCH 22/27] initial refinement: resize --- src/callbacks/refinement.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl index 03f7f6fe7..130de8756 100644 --- a/src/callbacks/refinement.jl +++ b/src/callbacks/refinement.jl @@ -41,6 +41,12 @@ function initial_refinement!(cb, 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 From 1e4959683e4a6de53f5cf5da5c1778df0f8eb3f6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:12:06 +0100 Subject: [PATCH 23/27] add split pattern --- src/refinement/refinement_pattern.jl | 93 ++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 src/refinement/refinement_pattern.jl diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl new file mode 100644 index 000000000..759fe8262 --- /dev/null +++ b/src/refinement/refinement_pattern.jl @@ -0,0 +1,93 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + 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) + end +end + +@inline function nchilds(system, refinement_pattern::CubicSplitting) + return 4 + refinement_pattern.center_particle +end + +struct TriangularSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + 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) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::TriangularSplitting) + return 3 + refinement_pattern.center_particle +end + +struct HexagonalSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + 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) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) + return 6 + refinement_pattern.center_particle +end From 278a83b114053e2c5b9810add7770adbcdcaabcf Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:19:32 +0100 Subject: [PATCH 24/27] add split method --- src/refinement/refinement.jl | 18 ++++-- src/refinement/split.jl | 103 +++++++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 6 deletions(-) create mode 100644 src/refinement/split.jl diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index ae435e481..35cf6de08 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,19 +1,24 @@ include("refinement_criteria.jl") -struct ParticleRefinement{RC, ELTYPE} - refinement_criteria :: RC - max_spacing_ratio :: ELTYPE - mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles +struct ParticleRefinement{SP, RC, ELTYPE} + splitting_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + n_particles_before_resize :: Int + n_new_particles :: Int end -function ParticleRefinement(; max_spacing_ratio, +function ParticleRefinement(; splitting_pattern, max_spacing_ratio, refinement_criteria=SpatialRefinementCriterion()) mass_ref = Vector{eltype(max_spacing_ratio)}() if refinement_criteria isa Tuple refinement_criteria = (refinement_criteria,) end - return ParticleRefinement(refinement_criteria, max_spacing_ratio, mass_ref) + + return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, + mass_ref, 0, 0) end resize_refinement!(system) = system @@ -37,6 +42,7 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) update_particle_spacing(semi, u_ode) # Split the particles (Algorithm 2) + split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) # Merge the particles (Algorithm 3) diff --git a/src/refinement/split.jl b/src/refinement/split.jl new file mode 100644 index 000000000..a04000b4e --- /dev/null +++ b/src/refinement/split.jl @@ -0,0 +1,103 @@ +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, refinement_pattern) = particle_refinement + + n_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 + n_split_candidates += 1 + end + end + + n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 + + particle_refinement.n_new_particles[] = n_split_candidates * n_childs_exclude_center + + 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 + (; mass_ref, max_spacing_ratio, refinement_pattern, n_particles_before_resize) = particle_refinement + (; alpha, relative_position) = refinement_pattern + + for particle in eachparticle(system) + m_a = hydrodynamic_mass(system, particle) + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a > m_max + smoothing_length_old = smoothing_length[particle] + mass_old = system.mass[particle] + + system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + + set_particle_pressure!(v, system, partice, pressure) + + set_particle_density!(v, system, partice, density) + + smoothing_length[particle] = alpha * smoothing_length_old + + pos_center = current_coords(system, u, particle) + vel_center = current_velocity(system, v, particle) + + for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) + child = n_particles_before_resize + child_id_local + + system.mass[child] = mass_old / nchilds(system, refinement_pattern) + + set_particle_pressure!(v, system, child, pressure) + + set_particle_density!(v, system, child, density) + + smoothing_length[child] = alpha * smoothing_length_old + + rel_pos = smoothing_length_old * relative_position[child] + new_pos = pos_center + rel_pos + + for dim in 1:ndims(system) + u[dim, child] = new_pos[dim] + v[dim, child] = vel_center[dim] + end + end + end + end + + return system +end From 94626984d9937a74ed45f98921833f95bbbee1a7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:27:56 +0100 Subject: [PATCH 25/27] include --- src/TrixiParticles.jl | 3 ++- src/refinement/refinement.jl | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d87a03796..54bfb7957 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 diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 35cf6de08..68d877d36 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,4 +1,6 @@ include("refinement_criteria.jl") +include("refinement_pattern.jl") +include("split.jl") struct ParticleRefinement{SP, RC, ELTYPE} splitting_pattern :: SP From 71c5347925a7ed7ae928e3dbabaf2a0a043ce69a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:23:30 +0100 Subject: [PATCH 26/27] rename to refinement_pattern --- src/refinement/refinement.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 68d877d36..0cf471cc8 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -3,7 +3,7 @@ include("refinement_pattern.jl") include("split.jl") struct ParticleRefinement{SP, RC, ELTYPE} - splitting_pattern :: SP + refinement_pattern :: SP refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles @@ -11,7 +11,7 @@ struct ParticleRefinement{SP, RC, ELTYPE} n_new_particles :: Int end -function ParticleRefinement(; splitting_pattern, max_spacing_ratio, +function ParticleRefinement(; refinement_pattern, max_spacing_ratio, refinement_criteria=SpatialRefinementCriterion()) mass_ref = Vector{eltype(max_spacing_ratio)}() @@ -19,7 +19,7 @@ function ParticleRefinement(; splitting_pattern, max_spacing_ratio, refinement_criteria = (refinement_criteria,) end - return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, + return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, mass_ref, 0, 0) end From 339678359f0b31db336693492970ee58c7bbfa44 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:33:38 +0100 Subject: [PATCH 27/27] fix `Ref` --- src/refinement/refinement.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 0cf471cc8..c81e53a0a 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -3,12 +3,12 @@ include("refinement_pattern.jl") include("split.jl") struct ParticleRefinement{SP, RC, ELTYPE} - refinement_pattern :: SP + refinement_pattern :: SP refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles - n_particles_before_resize :: Int - n_new_particles :: Int + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} end function ParticleRefinement(; refinement_pattern, max_spacing_ratio, @@ -20,7 +20,7 @@ function ParticleRefinement(; refinement_pattern, max_spacing_ratio, end return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, 0, 0) + mass_ref, Ref(0), Ref(0)) end resize_refinement!(system) = system