From 4ae0f8e6d1ad96d9d2a2e32d7e7b4d7464e1a376 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH] add h-factor --- .../fluid/entropically_damped_sph/rhs.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 96 ++++++++----------- src/schemes/fluid/fluid.jl | 22 +++++ .../fluid/weakly_compressible_sph/system.jl | 28 +----- 4 files changed, 69 insertions(+), 79 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 7b230ddc3..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -63,53 +63,53 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, 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)), - 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) - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - mass = copy(initial_condition.mass) - - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end + 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) + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + mass = copy(initial_condition.mass) + + if ndims(smoothing_kernel) != NDIMS + throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) + end - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - nothing) + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, + density_calculator, + NDIMS, ELTYPE, + nothing) - nu_edac = (alpha * smoothing_length * sound_speed) / 8 + nu_edac = (alpha * smoothing_length * sound_speed) / 8 - cache = create_cache_density(initial_condition, density_calculator) - cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) - cache = (; - create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, - nparticles(initial_condition))..., - cache...) + cache = create_cache_density(initial_condition, density_calculator) + cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) - return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, smoothing_kernel, - sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, transport_velocity, source_terms, - buffer, particle_refinement, cache) - end + return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, + smoothing_kernel, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, + transport_velocity, source_terms, buffer, + particle_refinement, cache) +end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) @nospecialize system # reduce precompilation time @@ -147,18 +147,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