Adaptive Particle Refinement: Variable smoothing length

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

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

src/general/semidiscretization.jl

@@ -134,8 +134,8 @@ function create_neighborhood_search(neighborhood_search, system, neighbor)
 end
 
 @inline function compact_support(system, neighbor)
-    (; smoothing_kernel, smoothing_length) = system
-    return compact_support(smoothing_kernel, smoothing_length)
+    (; smoothing_kernel) = system
+    return compact_support(smoothing_kernel, smoothing_length(system, 1)) # TODO
 end
 
 @inline function compact_support(system::OpenBoundarySPHSystem, neighbor)
@@ -407,7 +407,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)

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.

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
                            +

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)

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

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -44,12 +44,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                     gravity-like source terms.
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
-                                   PF, ST, B, C} <: FluidSystem{NDIMS, IC}
+                                   PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC}
     initial_condition                 :: IC
     mass                              :: M # Vector{ELTYPE}: [particle]
     density_calculator                :: DC
     smoothing_kernel                  :: K
-    smoothing_length                  :: ELTYPE
     sound_speed                       :: ELTYPE
     viscosity                         :: V
     nu_edac                           :: ELTYPE
@@ -59,55 +58,57 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
     transport_velocity                :: TV
     source_terms                      :: ST
     buffer                            :: B
+    particle_refinement               :: PR
     cache                             :: C
+end
 
-    function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
-                                         smoothing_length, sound_speed;
-                                         pressure_acceleration=inter_particle_averaged_pressure,
-                                         density_calculator=SummationDensity(),
-                                         transport_velocity=nothing,
-                                         alpha=0.5, viscosity=nothing,
-                                         acceleration=ntuple(_ -> 0.0,
-                                                             ndims(smoothing_kernel)),
-                                         source_terms=nothing, buffer_size=nothing)
-        buffer = isnothing(buffer_size) ? nothing :
-                 SystemBuffer(nparticles(initial_condition), buffer_size)
+function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
+                                     smoothing_length, sound_speed;
+                                     pressure_acceleration=inter_particle_averaged_pressure,
+                                     density_calculator=SummationDensity(),
+                                     transport_velocity=nothing,
+                                     alpha=0.5, viscosity=nothing,
+                                     acceleration=ntuple(_ -> 0.0,
+                                                         ndims(smoothing_kernel)),
+                                     particle_refinement=nothing,
+                                     source_terms=nothing, buffer_size=nothing)
+    buffer = isnothing(buffer_size) ? nothing :
+             SystemBuffer(nparticles(initial_condition), buffer_size)
 
-        initial_condition = allocate_buffer(initial_condition, buffer)
+    initial_condition = allocate_buffer(initial_condition, buffer)
 
-        NDIMS = ndims(initial_condition)
-        ELTYPE = eltype(initial_condition)
+    NDIMS = ndims(initial_condition)
+    ELTYPE = eltype(initial_condition)
 
-        mass = copy(initial_condition.mass)
+    mass = copy(initial_condition.mass)
 
-        if ndims(smoothing_kernel) != NDIMS
-            throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
-        end
+    if ndims(smoothing_kernel) != NDIMS
+        throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
+    end
 
-        acceleration_ = SVector(acceleration...)
-        if length(acceleration_) != NDIMS
-            throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
-        end
+    acceleration_ = SVector(acceleration...)
+    if length(acceleration_) != NDIMS
+        throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
+    end
 
-        pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
-                                                                         density_calculator,
-                                                                         NDIMS, ELTYPE,
-                                                                         nothing)
+    pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
+                                                                     density_calculator,
+                                                                     NDIMS, ELTYPE,
+                                                                     nothing)
 
-        nu_edac = (alpha * smoothing_length * sound_speed) / 8
+    nu_edac = (alpha * smoothing_length * sound_speed) / 8
 
-        cache = create_cache_density(initial_condition, density_calculator)
-        cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
+    cache = create_cache_density(initial_condition, density_calculator)
+    cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
+    cache = (;
+             create_cache_refinement(initial_condition, particle_refinement,
+                                     smoothing_length)..., cache...)
 
-        new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass),
-            typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity),
-            typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms),
-            typeof(buffer),
-            typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel,
-                           smoothing_length, sound_speed, viscosity, nu_edac, acceleration_,
-                           nothing, pressure_acceleration, transport_velocity, source_terms,
-                           buffer, cache)
-    end
+    return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator,
+                                       smoothing_kernel, sound_speed, viscosity, nu_edac,
+                                       acceleration_, nothing, pressure_acceleration,
+                                       transport_velocity, source_terms, buffer,
+                                       particle_refinement, cache)
 end
 
 function Base.show(io::IO, system::EntropicallyDampedSPHSystem)

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
 

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

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