diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 4689920c3..1d1450d8b 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -203,18 +203,16 @@ function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing return zero(initial_density) end -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, +@inline function pressure_acceleration(pressure_correction, m_b, + particle, boundary_particle, + particle_system, boundary_system, boundary_model::BoundaryModelDummyParticles, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator) (; density_calculator) = boundary_model - pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, v_boundary_system, + pressure_acceleration(pressure_correction, m_b, particle, boundary_particle, + particle_system, boundary_system, boundary_model, density_calculator, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator) @@ -223,17 +221,13 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with ContinuityDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, +@inline function pressure_acceleration(pressure_correction, m_b, + particle, boundary_particle, + particle_system, boundary_system, boundary_model::BoundaryModelDummyParticles, boundary_density_calculator, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator::ContinuityDensity) - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) - return -m_b * (particle_system.pressure[particle] + boundary_model.pressure[boundary_particle]) / (rho_a * rho_b) * grad_kernel @@ -242,17 +236,13 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with SummationDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, +@inline function pressure_acceleration(pressure_correction, m_b, + particle, boundary_particle, + particle_system, boundary_system, boundary_model::BoundaryModelDummyParticles, boundary_density_calculator, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator::SummationDensity) - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) - return -m_b * (particle_system.pressure[particle] / rho_a^2 + boundary_model.pressure[boundary_particle] / rho_b^2) * @@ -262,17 +252,13 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with ContinuityDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, +@inline function pressure_acceleration(pressure_correction, m_b, + particle, boundary_particle, + particle_system, boundary_system, boundary_model::BoundaryModelDummyParticles, ::PressureMirroring, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator::ContinuityDensity) - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) - return -m_b * (particle_system.pressure[particle] + particle_system.pressure[particle]) / (rho_a * rho_b) * grad_kernel @@ -281,17 +267,13 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with SummationDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, +@inline function pressure_acceleration(pressure_correction, m_b, + particle, boundary_particle, + particle_system, boundary_system, boundary_model::BoundaryModelDummyParticles, ::PressureMirroring, rho_a, rho_b, pos_diff, distance, grad_kernel, fluid_density_calculator::SummationDensity) - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) - return -m_b * (particle_system.pressure[particle] / rho_a^2 + particle_system.pressure[particle] / rho_b^2) * diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index bda187b82..758d5559c 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -76,9 +76,8 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, neighbor_system, - v_neighbor_system, +@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system, boundary_model::BoundaryModelMonaghanKajtar, rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 9d7c67921..d5c47e605 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,15 +6,15 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) - (; density_calculator, state_equation, correction, smoothing_length) = particle_system + (; density_calculator, state_equation, correction) = particle_system (; sound_speed) = state_equation viscosity = viscosity_model(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code - # and the two other lines below that are marked as "debug example". + # In order to visualize quantities like pressure forces or viscosity forces, uncomment + # the following code and the two other lines below that are marked as "debug example". # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) # Loop over all pairs of particles and neighbors within the kernel cutoff. @@ -38,10 +38,9 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, - particle_system, v_particle_system, - neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, pos_diff, + dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system, + rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) dv_viscosity = viscosity_correction * @@ -74,10 +73,10 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with ContinuityDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, +@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system::WeaklyCompressibleSPHSystem, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, + rho_a, rho_b, pos_diff, distance, grad_kernel, ::ContinuityDensity) return (-m_b * (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) / @@ -88,10 +87,10 @@ end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with SummationDensity. -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, +@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system::WeaklyCompressibleSPHSystem, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, + rho_a, rho_b, pos_diff, distance, grad_kernel, ::SummationDensity) return (-m_b * (particle_system.pressure[particle] / rho_a^2 + @@ -99,17 +98,17 @@ end pressure_correction end -@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, +@inline function pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, + rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) (; boundary_model) = neighbor_system - return pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, neighbor_system, - v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, + return pressure_acceleration(pressure_correction, m_b, particle, neighbor, + particle_system, neighbor_system, + boundary_model, rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 163388519..ccf422a44 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -98,9 +98,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Boundary forces # Note: neighbor and particle are switched in this call and `pressure_correction` is set to `1.0` (no correction) - dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system, - v_neighbor_system, particle, - particle_system, v_particle_system, + dv_boundary = pressure_acceleration(1.0, m_b, neighbor, particle, + neighbor_system, particle_system, boundary_model, rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index bf1d4b5eb..e7ee22392 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -1,5 +1,3 @@ -using QuadGK - @testset verbose=true "Smoothing Kernels" begin # Don't show all kernel tests in the final overview @testset verbose=false "Integral" begin diff --git a/test/schemes/fluid/weakly_compressible_sph/rhs.jl b/test/schemes/fluid/weakly_compressible_sph/rhs.jl new file mode 100644 index 000000000..f9c441fea --- /dev/null +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -0,0 +1,166 @@ +@testset verbose=true "WCSPH RHS" begin + @testset verbose=true "`pressure_acceleration`" begin + # Use `@trixi_testset` to isolate the mock functions in a separate namespace + @trixi_testset "Symmetry" begin + density_calculators = [ContinuityDensity(), SummationDensity()] + masses = [[0.01, 0.01], [0.73, 0.31]] + densities = [ + [1000.0, 1000.0], + [1000.0, 1000.0], + [900.0, 1201.0], + [1003.0, 353.4], + ] + pressures = [ + [0.0, 0.0], + [10_000.0, 10_000.0], + [10.0, 10_000.0], + [1000.0, -1000.0], + ] + grad_kernels = [0.3, 104.0] + particle = 2 + neighbor = 3 + + # Not used for fluid-fluid interaction + pos_diff = 0 + distance = 0 + + @testset "`$(nameof(typeof(density_calculator)))`" for density_calculator in density_calculators + for (m_a, m_b) in masses, (rho_a, rho_b) in densities, + (p_a, p_b) in pressures, grad_kernel in grad_kernels + + # Partly copied from constructor test, just to create a WCSPH system + coordinates = zeros(2, 3) + velocity = zeros(2, 3) + mass = zeros(3) + density = zeros(3) + state_equation = Val(:state_equation) + smoothing_kernel = Val(:smoothing_kernel) + TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 + smoothing_length = -1.0 + + initial_condition = InitialCondition(coordinates, velocity, mass, + density) + system = WeaklyCompressibleSPHSystem(initial_condition, + density_calculator, + state_equation, smoothing_kernel, + smoothing_length) + + # `system` is only used for the pressure + system.pressure .= [0.0, p_a, p_b] + + # Compute accelerations a -> b and b -> a + dv1 = TrixiParticles.pressure_acceleration(1.0, m_b, particle, neighbor, + system, system, rho_a, rho_b, + pos_diff, distance, + grad_kernel, + density_calculator) + + dv2 = TrixiParticles.pressure_acceleration(1.0, m_a, neighbor, particle, + system, system, rho_b, rho_a, + -pos_diff, distance, + -grad_kernel, + density_calculator) + + # Test that both forces are identical but in opposite directions + @test isapprox(m_a * dv1, -m_b * dv2, rtol=2eps()) + end + end + end + end + + @testset verbose=true "`interact!`" begin + # The following tests for linear and angular momentum conservation + # are based on Section 3.3.4 of + # Daniel J. Price. "Smoothed Particle Hydrodynamics and Magnetohydrodynamics." + # In: Journal of Computational Physics 231.3 (2012), pages 759–94. + # https://doi.org/10.1016/j.jcp.2010.12.011 + @testset verbose=true "Momentum Conservation" begin + # We are testing the momentum conservation of SPH with random initial configurations + density_calculators = [ContinuityDensity(), SummationDensity()] + + # Random initial configuration + mass = [ + [3.11, 1.55, 2.22, 3.48, 0.21, 3.73, 0.21, 3.45], + [0.82, 1.64, 1.91, 0.02, 0.08, 1.58, 4.94, 0.7], + ] + density = [ + [914.34, 398.36, 710.22, 252.54, 843.81, 694.73, 670.5, 539.14], + [280.15, 172.25, 267.1, 130.42, 382.3, 477.21, 848.31, 188.62], + ] + pressure = [ + [91438.0, 16984.0, 58638.0, 10590.0, 92087.0, 66586.0, 64723.0, 49862.0], + [31652.0, -21956.0, 2874.0, -12489.0, 27206.0, 32225.0, 42848.0, 3001.0], + ] + coordinates = [ + [0.16 0.55 0.08 0.58 0.52 0.26 0.32 0.99; + 0.76 0.6 0.47 0.4 0.25 0.79 0.45 0.63], + [0.4 0.84 0.47 0.02 0.64 0.85 0.02 0.15; + 0.83 0.62 0.99 0.57 0.25 0.72 0.34 0.69], + ] + velocity = [ + [1.05 0.72 0.12 1.22 0.67 0.85 1.42 -0.57; + 1.08 0.68 0.74 -0.27 -1.22 0.43 1.41 1.25], + [-1.84 -1.4 5.21 -5.99 -5.02 9.5 -4.51 -8.28; + 0.78 0.1 9.67 8.46 9.29 5.18 -4.83 -4.87]] + + # The state equation is only needed to unpack `sound_speed`, so we can mock + # it by using a `NamedTuple`. + state_equation = (; sound_speed=0.0) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.3 + search_radius = TrixiParticles.compact_support(smoothing_kernel, + smoothing_length) + + @testset "`$(nameof(typeof(density_calculator)))`" for density_calculator in density_calculators + for i in eachindex(mass) + initial_condition = InitialCondition(coordinates[i], velocity[i], + mass[i], density[i]) + system = WeaklyCompressibleSPHSystem(initial_condition, + density_calculator, + state_equation, smoothing_kernel, + smoothing_length) + + # Overwrite `system.pressure` + system.pressure .= pressure[i] + + u = coordinates[i] + if density_calculator isa SummationDensity + # Density is stored in the cache + v = velocity[i] + system.cache.density .= density[i] + else + # Density is integrated with `ContinuityDensity` + v = vcat(velocity[i], density[i]') + end + + nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius, + TrixiParticles.eachparticle(system)) + + # Result + dv = zeros(3, 8) + TrixiParticles.interact!(dv, v, u, v, u, nhs, system, system) + + # Linear momentum conservation + # ∑ m_a dv_a + deriv_linear_momentum = sum(mass[i]' .* view(dv, 1:2, :), dims=2) + + @test isapprox(deriv_linear_momentum, zeros(2, 1), atol=3e-14) + + # Angular momentum conservation + # m_a (r_a × dv_a) + function deriv_angular_momentum(particle) + r_a = SVector(u[1, particle], u[2, particle], 0.0) + dv_a = SVector(dv[1, particle], dv[2, particle], 0.0) + + return mass[i][particle] * cross(r_a, dv_a) + end + + # ∑ m_a (r_a × dv_a) + deriv_angular_momentum = sum(deriv_angular_momentum, 1:8) + + @test isapprox(deriv_angular_momentum, zeros(3), atol=2e-14) + end + end + end + end +end diff --git a/test/schemes/fluid/weakly_compressible_sph/weakly_compressible_sph.jl b/test/schemes/fluid/weakly_compressible_sph/weakly_compressible_sph.jl index 466889a5c..110143ae1 100644 --- a/test/schemes/fluid/weakly_compressible_sph/weakly_compressible_sph.jl +++ b/test/schemes/fluid/weakly_compressible_sph/weakly_compressible_sph.jl @@ -1,2 +1,3 @@ include("density_diffusion.jl") include("state_equation.jl") +include("rhs.jl") diff --git a/test/test_util.jl b/test/test_util.jl index 9d714b727..6897bad29 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -2,6 +2,7 @@ using Test using TrixiParticles using LinearAlgebra using Printf +using QuadGK: quadgk """ @trixi_testset "name of the testset" #= code to test #=