add h-factor

trixi-framework · Nov 18, 2024 · 4ae0f8e · 4ae0f8e
1 parent e99571f
commit 4ae0f8e
Show file tree

Hide file tree

Showing 4 changed files with 69 additions and 79 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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