From afb6e9ffed0559fe645820242207a26e02299327 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 1/3] 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 cf6c8bed6..a82bdf8ae 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 95c9ff88e..9af1f9ebc 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 5d712235f..d6dac2f36 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()) && 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 e99571f21bbb3bf7f85d2978ba3b21f38f18098d 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 2/3] Make it work with EDAC --- src/general/semidiscretization.jl | 8 ++-- .../dummy_particles/dummy_particles.jl | 6 ++- .../fluid/entropically_damped_sph/system.jl | 39 ++++++++++++------- .../fluid/weakly_compressible_sph/system.jl | 2 +- src/visualization/write2vtk.jl | 2 +- 5 files changed, 38 insertions(+), 19 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..6c0c7e37b 100644 --- a/src/general/semidiscretization.jl +++ b/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) 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 eb2aa1c1d..7b230ddc3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -44,12 +44,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K - smoothing_length :: ELTYPE sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -59,9 +58,11 @@ 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, +function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), @@ -69,6 +70,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, 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) @@ -98,17 +100,16 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, 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) + cache = (; + create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, + nparticles(initial_condition))..., + 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 -end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) @nospecialize system # reduce precompilation time @@ -146,6 +147,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 4ae0f8e6d1ad96d9d2a2e32d7e7b4d7464e1a376 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 3/3] 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