From 87b0ea124cec625bbb3322237651bc70325fad4d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 Dec 2023 16:42:44 +0100 Subject: [PATCH 01/12] Fix Issue #245 (#314) * implement * Revert "implement" This reverts commit 171505a4e078c7271e838c35f98f3df7f4cf77de. * fix Check if PVD file already contains older data #245 * cleanup * format * fix --- src/visualization/write2vtk.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index c6eaace56..d23f0b9a9 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -66,7 +66,8 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix collection_file = joinpath(output_directory, add_opt_str_pre(prefix) * "$system_name") - pvd = paraview_collection(collection_file; append=true) + # Reset the collection when the iteration is 0 + pvd = paraview_collection(collection_file; append=iter > 0) points = periodic_coords(current_coordinates(u, system), periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] From 39141b0572e4436a152732666097cc3573ee365f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 21 Dec 2023 10:27:12 +0100 Subject: [PATCH 02/12] Typo (#323) Rectangular Tank does not work as expected --- src/setups/rectangular_tank.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index e38f5cdee..cafd439a9 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -152,7 +152,7 @@ function fluid_particles_per_dimension(size::NTuple{3}, particle_spacing) "fluid length in x-direction") n_particles_y, new_y_size = round_n_particles(size[2], particle_spacing, "fluid length in y-direction") - n_particles_z, new_z_size = round_n_particles(size[2], particle_spacing, + n_particles_z, new_z_size = round_n_particles(size[3], particle_spacing, "fluid length in z-direction") return (n_particles_x, n_particles_y, n_particles_z), From f43f15621d690ead4089a77eeb4497b108dc0403 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Dec 2023 16:54:41 +0100 Subject: [PATCH 03/12] Another fix for RISC processors (#325) * fix * format --- .../boundary/dummy_particles/dummy_particles.jl | 14 ++++++++++---- .../fluid/weakly_compressible_sph/system.jl | 3 ++- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index c967b827c..a99387cdf 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -319,18 +319,24 @@ end function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityDensity}, system, v, u, v_ode, u_ode, semi) - (; state_equation, pressure) = boundary_model # Limit pressure to be non-negative to avoid attractive forces between fluid and # boundary particles at free surfaces (sticking artifacts). - @trixi_timeit timer() "state equation" @threaded for particle in eachparticle(system) - pressure[particle] = max(state_equation(particle_density(v, boundary_model, - particle)), 0.0) + @threaded for particle in eachparticle(system) + apply_state_equation!(boundary_model, particle_density(v, boundary_model, + particle), particle) end return boundary_model end +# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). +# Otherwise, `@threaded` does not work here with Julia ARM on macOS. +# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. +@inline function apply_state_equation!(boundary_model, density, particle) + boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) +end + function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, v_ode, u_ode, semi) (; pressure, state_equation, cache, viscosity) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 778999913..e65d1465a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -251,7 +251,8 @@ end # Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). # Otherwise, `@threaded` does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. -@inline function apply_state_equation!(system, density, particle) +@inline function apply_state_equation!(system::WeaklyCompressibleSPHSystem, density, + particle) system.pressure[particle] = system.state_equation(density) end From 2a288da7d512741c3042c62379826992419c95d8 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Dec 2023 16:54:59 +0100 Subject: [PATCH 04/12] Another fix for RISC processors (#325) * fix * format From 205d69667f2f8161b89809028885d7f2f8ac3af7 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 2 Jan 2024 11:36:21 +0100 Subject: [PATCH 05/12] Add common files which could be accidentally pushed (#320) * add common files which could be accidentally pushed * [skip ci] --- .gitignore | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.gitignore b/.gitignore index b83d9967c..e017a301c 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,10 @@ docs/src/license.md docs/src/code_of_conduct.md docs/src/news.md run +*.json +*.png +*.gif +*.svg +*.vtu +*.pvd +*.mp4 From 4101f1169fd6bb2a7069c0ae573a0405ba11b863 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Tue, 2 Jan 2024 15:01:23 +0100 Subject: [PATCH 06/12] `min_coordinates` for `RectangularTank` (#298) * add `min_coordinates` * add tests * indentation * implement suggestions --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/general/semidiscretization.jl | 9 ++--- src/setups/rectangular_tank.jl | 7 ++++ test/setups/rectangular_tank.jl | 57 +++++++++++++++++++------------ 3 files changed, 45 insertions(+), 28 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6c3bf4796..6cfe1f8d2 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -168,8 +168,7 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`. function semidiscretize(semi, tspan; reset_threads=true) (; systems) = semi - @assert all(system -> eltype(system) === eltype(systems[1]), - systems) + @assert all(system -> eltype(system) === eltype(systems[1]), systems) ELTYPE = eltype(systems[1]) # Optionally reset Polyester.jl threads. See @@ -190,10 +189,8 @@ function semidiscretize(semi, tspan; reset_threads=true) end end - sizes_u = (u_nvariables(system) * n_moving_particles(system) - for system in systems) - sizes_v = (v_nvariables(system) * n_moving_particles(system) - for system in systems) + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in systems) + sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in systems) u0_ode = Vector{ELTYPE}(undef, sum(sizes_u)) v0_ode = Vector{ELTYPE}(undef, sum(sizes_v)) diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index cafd439a9..aaf013ada 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -1,6 +1,7 @@ @doc raw""" RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density; n_layers=1, spacing_ratio=1.0, + min_coordinates=zeros(length(fluid_size)), init_velocity=zeros(length(fluid_size)), boundary_density=fluid_density, faces=Tuple(trues(2 * length(fluid_size)))) @@ -18,6 +19,7 @@ Rectangular tank filled with a fluid to set up dam-break-style simulations. - `spacing_ratio`: Ratio of `particle_spacing` to boundary particle spacing. A value of 2 means that the boundary particle spacing will be half the fluid particle spacing. +- `min_coordinates`: Coordinates of the corner in negative coordinate directions. - `init_velocity`: The initial velocity of each fluid particle as `(x, y)` (or `(x, y, z)` in 3D). - `boundary_density`: Density of each boundary particle (by default set to the rest density) - `faces`: By default all faces are generated. Set faces by passing a @@ -58,6 +60,7 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} function RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density; pressure=0.0, n_layers=1, spacing_ratio=1.0, + min_coordinates=zeros(length(fluid_size)), init_velocity=zeros(length(fluid_size)), boundary_density=fluid_density, faces=Tuple(trues(2 * length(fluid_size))), @@ -119,6 +122,10 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} boundary_masses, boundary_densities, particle_spacing=boundary_spacing) + # Move the tank corner in the negative coordinate directions to the desired position + fluid.coordinates .+= min_coordinates + boundary.coordinates .+= min_coordinates + return new{NDIMS, 2 * NDIMS, ELTYPE}(fluid, boundary, fluid_size_, tank_size_, faces, face_indices, particle_spacing, spacing_ratio, n_layers, diff --git a/test/setups/rectangular_tank.jl b/test/setups/rectangular_tank.jl index 278cc6951..8b9932101 100644 --- a/test/setups/rectangular_tank.jl +++ b/test/setups/rectangular_tank.jl @@ -10,6 +10,7 @@ @testset "Coordinates" begin particle_spacings = [0.1, 0.2] spacing_ratios = [1, 3] + min_coordinates = [(0.0, 0.0), (-0.3, 2.0)] expected_fluid_coords = [ [0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45; @@ -33,17 +34,22 @@ ], ] - @testset "Particle Spacing: $(particle_spacings[i])" for i in eachindex(particle_spacings) - @testset "Spacing Ratio: $(spacing_ratios[j])" for j in eachindex(spacing_ratios) - tank = RectangularTank(particle_spacings[i], - (water_width, water_height), - (tank_width, tank_height), - water_density, - spacing_ratio=spacing_ratios[j]) - - @test isapprox(tank.fluid.coordinates, expected_fluid_coords[i]) - @test isapprox(tank.boundary.coordinates, - expected_bound_coords[i][j]) + @testset "Move Tank: $(min_coordinates[i])" for i in eachindex(min_coordinates) + @testset "Particle Spacing: $(particle_spacings[j])" for j in eachindex(particle_spacings) + @testset "Spacing Ratio: $(spacing_ratios[k])" for k in eachindex(spacing_ratios) + tank = RectangularTank(particle_spacings[j], + (water_width, water_height), + (tank_width, tank_height), + water_density, + spacing_ratio=spacing_ratios[k], + min_coordinates=min_coordinates[i]) + expected_fluid_coords_ = copy(expected_fluid_coords[j]) + expected_bound_coords_ = copy(expected_bound_coords[j][k]) + expected_fluid_coords_ .+= min_coordinates[i] + expected_bound_coords_ .+= min_coordinates[i] + @test isapprox(tank.fluid.coordinates, expected_fluid_coords_) + @test isapprox(tank.boundary.coordinates, expected_bound_coords_) + end end end end @@ -267,6 +273,7 @@ end @testset "Coordinates" begin particle_spacings = [0.1, 0.2] spacing_ratios = [1, 3] + min_coordinates = [(0.0, 0.0, 0.0), (-0.3, 2.0, -0.5)] expected_fluid_coords = [ [0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35; @@ -296,17 +303,23 @@ end ], ] - @testset "Particle Spacing: $(particle_spacings[i])" for i in eachindex(particle_spacings) - @testset "Spacing Ratio: $(spacing_ratios[j])" for j in eachindex(spacing_ratios) - tank = RectangularTank(particle_spacings[i], - (water_width, water_height, water_depth), - (tank_width, tank_height, tank_depth), - water_density, - spacing_ratio=spacing_ratios[j]) - - @test isapprox(tank.fluid.coordinates, expected_fluid_coords[i]) - @test isapprox(tank.boundary.coordinates, - expected_bound_coords[i][j]) + @testset "Move Tank: $(min_coordinates[i])" for i in eachindex(min_coordinates) + @testset "Particle Spacing: $(particle_spacings[j])" for j in eachindex(particle_spacings) + @testset "Spacing Ratio: $(spacing_ratios[k])" for k in eachindex(spacing_ratios) + tank = RectangularTank(particle_spacings[j], + (water_width, water_height, water_depth), + (tank_width, tank_height, tank_depth), + water_density, + spacing_ratio=spacing_ratios[k], + min_coordinates=min_coordinates[i]) + expected_fluid_coords_ = copy(expected_fluid_coords[j]) + expected_bound_coords_ = copy(expected_bound_coords[j][k]) + expected_fluid_coords_ .+= min_coordinates[i] + expected_bound_coords_ .+= min_coordinates[i] + + @test isapprox(tank.fluid.coordinates, expected_fluid_coords_) + @test isapprox(tank.boundary.coordinates, expected_bound_coords_) + end end end end From 4876752913fb743700ce2c4dafc84e2e34c473b7 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Wed, 3 Jan 2024 10:38:49 +0100 Subject: [PATCH 07/12] Circle segment support for `RoundSphere` (#297) * add `theta` * add angular range and docs * add tests * fix typo * remove utf symbol * implement suggestions --- src/setups/sphere_shape.jl | 52 +++++++++++++++++++++++++++++-------- test/setups/sphere_shape.jl | 6 +++++ 2 files changed, 47 insertions(+), 11 deletions(-) diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index 73a40a376..a2ab46b9c 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -127,15 +127,38 @@ as it will have corners (like a sphere in Minecraft). struct VoxelSphere end """ - RoundSphere() + RoundSphere(; start_angle=0.0, end_angle=2π) -Construct a sphere by nesting perfectly round concentric spheres. +Construct a sphere (or sphere segment) by nesting perfectly round concentric spheres. The resulting ball will be perfectly round, but will not have a regular inner structure. +# Keywords +- `start_angle`: The starting angle of the sphere segment in radians. It determines the + beginning point of the segment. The default is set to `0.0` representing + the positive x-axis. +- `end_angle`: The ending angle of the sphere segment in radians. It defines the termination + point of the segment. The default is set to `2pi`, completing a full sphere. + !!! note "Usage" See [`SphereShape`](@ref) on how to use this. + +!!! warning "Warning" + The sphere segment is intended for 2D geometries and hollow spheres. If used for filled + spheres or in a 3D context, results may not be accurate. """ -struct RoundSphere end +struct RoundSphere{AR} + angle_range::AR + + function RoundSphere(; start_angle=0.0, end_angle=2pi) + if start_angle > end_angle + throw(ArgumentError("`end_angle` should be greater than `start_angle`")) + end + + angle_range = (start_angle, end_angle) + + new{typeof(angle_range)}(angle_range) + end +end function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_position, n_layers, layer_outwards, tlsph) @@ -193,7 +216,7 @@ function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_pos return reinterpret(reshape, ELTYPE, coords) end -function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, +function sphere_shape_coords(sphere::RoundSphere, particle_spacing, radius, center, n_layers, layer_outwards, tlsph) if n_layers > 0 if layer_outwards @@ -227,7 +250,7 @@ function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, coords = zeros(length(center), 0) for layer in 0:(n_layers - 1) - sphere_coords = round_sphere(particle_spacing, + sphere_coords = round_sphere(sphere, particle_spacing, inner_radius + layer * particle_spacing, center) coords = hcat(coords, sphere_coords) end @@ -235,8 +258,11 @@ function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, return coords end -function round_sphere(particle_spacing, radius, center::SVector{2}) - n_particles = round(Int, 2pi * radius / particle_spacing) +function round_sphere(sphere, particle_spacing, radius, center::SVector{2}) + (; angle_range) = sphere + + theta = angle_range[2] - angle_range[1] + n_particles = round(Int, theta * radius / particle_spacing) if n_particles <= 2 # 2 or less particles produce weird, asymmetric results. @@ -244,8 +270,12 @@ function round_sphere(particle_spacing, radius, center::SVector{2}) return collect(reshape(center, (2, 1))) end - # Remove the last particle at 2pi, which overlaps with the first at 0 - t = LinRange(0, 2pi, n_particles + 1)[1:(end - 1)] + if !isapprox(theta, 2pi) + t = LinRange(angle_range[1], angle_range[2], n_particles + 1) + else + # Remove the last particle at 2pi, which overlaps with the first at 0 + t = LinRange(angle_range[1], angle_range[2], n_particles + 1)[1:(end - 1)] + end particle_coords = Array{Float64, 2}(undef, 2, length(t)) @@ -256,7 +286,7 @@ function round_sphere(particle_spacing, radius, center::SVector{2}) return particle_coords end -function round_sphere(particle_spacing, radius, center::SVector{3}) +function round_sphere(sphere, particle_spacing, radius, center::SVector{3}) # The number of particles can either be calculated in 2D or in 3D. # Let δ be the particle spacing and r the sphere radius. # @@ -374,7 +404,7 @@ function round_sphere(particle_spacing, radius, center::SVector{3}) circle_spacing = 1.0 end - circle_coords_2d = round_sphere(circle_spacing, circle_radius, + circle_coords_2d = round_sphere(sphere, circle_spacing, circle_radius, SVector(center[1], center[2])) circle_coords_3d = vcat(circle_coords_2d, center[3] .+ z * ones(1, size(circle_coords_2d, 2))) diff --git a/test/setups/sphere_shape.jl b/test/setups/sphere_shape.jl index e39cb2151..7577ed1ab 100644 --- a/test/setups/sphere_shape.jl +++ b/test/setups/sphere_shape.jl @@ -16,6 +16,7 @@ "Hollow RoundSphere with One Layer", "Hollow RoundSphere with Multiple Layers", "Hollow RoundSphere Outwards", + "Hollow RoundSphere Segment", ] shapes = [ @@ -39,6 +40,9 @@ SphereShape(particle_spacing, radius, (-3.0, 4.0), 1000.0, sphere_type=RoundSphere(), n_layers=3, layer_outwards=true), + SphereShape(particle_spacing, radius, (3.0, -4), 1000.0, + sphere_type=RoundSphere(start_angle=pi, end_angle=1.5pi), + n_layers=3, layer_outwards=true), ] expected_coords = [ @@ -66,6 +70,8 @@ -4.0 -3.22865487 -2.8182307 -2.96076952 -3.58957583 -4.41042417 -5.03923048 -5.1817693 -4.77134513 -4.0 -3.23463314 -2.58578644 -2.15224093 -2.0 -2.15224093 -2.58578644 -3.23463314 -4.0 -4.76536686 -5.41421356 -5.84775907 -6.0 -5.84775907 -5.41421356 -4.76536686], [-0.2 -0.31341967 -0.64449011 -1.16638994 -1.83683796 -2.60151845 -3.39848155 -4.16316204 -4.83361006 -5.35550989 -5.68658033 -5.8 -5.68658033 -5.35550989 -4.83361006 -4.16316204 -3.39848155 -2.60151845 -1.83683796 -1.16638994 -0.64449011 -0.31341967 0.6 0.50974048 0.24348792 -0.18540666 -0.75543671 -1.43801854 -2.19892464 -3.0 -3.80107536 -4.56198146 -5.24456329 -5.81459334 -6.24348792 -6.50974048 -6.6 -6.50974048 -6.24348792 -5.81459334 -5.24456329 -4.56198146 -3.80107536 -3.0 -2.19892464 -1.43801854 -0.75543671 -0.18540666 0.24348792 0.50974048 1.4 1.32929019 1.11943343 0.77717469 0.31351445 -0.25664487 -0.91497789 -1.64032522 -2.40937363 -3.19740525 -3.97909211 -4.72931014 -5.42394672 -6.04067566 -6.55967478 -6.96426302 -7.24143659 -7.38228689 -7.38228689 -7.24143659 -6.96426302 -6.55967478 -6.04067566 -5.42394672 -4.72931014 -3.97909211 -3.19740525 -2.40937363 -1.64032522 -0.91497789 -0.25664487 0.31351445 0.77717469 1.11943343 1.32929019; 4.0 4.78885116 5.51379429 6.11609881 6.54696959 6.77150004 6.77150004 6.54696959 6.11609881 5.51379429 4.78885116 4.0 3.21114884 2.48620571 1.88390119 1.45303041 1.22849996 1.22849996 1.45303041 1.88390119 2.48620571 3.21114884 4.0 4.80107536 5.56198146 6.24456329 6.81459334 7.24348792 7.50974048 7.6 7.50974048 7.24348792 6.81459334 6.24456329 5.56198146 4.80107536 4.0 3.19892464 2.43801854 1.75543671 1.18540666 0.75651208 0.49025952 0.4 0.49025952 0.75651208 1.18540666 1.75543671 2.43801854 3.19892464 4.0 4.78565034 5.54604923 6.25675682 6.89493039 7.44005852 7.87462034 8.18464867 8.36017895 8.39556949 8.28968281 8.0459222 7.67212232 7.1802974 6.58625511 5.90908845 5.17056212 4.39441296 3.60558704 2.82943788 2.09091155 1.41374489 0.8197026 0.32787768 -0.0459222 -0.28968281 -0.39556949 -0.36017895 -0.18464867 0.12537966 0.55994148 1.10506961 1.74324318 2.45395077 3.21434966], + [0.20000000000000018 0.33704175437357 0.7347524157501475 1.354201293581075 2.134752415750147 2.9999999999999996 -0.5999999999999996 -0.5097404838545647 -0.2434879244487087 0.18540666311509257 0.7554367133085589 1.4380185391767903 2.1989246377572673 2.9999999999999996 -1.4000000000000004 -1.3331541132537161 -1.1346475314579978 -0.8105117766515302 -0.37059554972350384 0.17173451737922596 0.8000000000000016 1.4951113693670546 2.2359480182655065 2.999999999999999; + -3.9999999999999996 -4.865247584249852 -5.645798706418924 -6.2652475842498525 -6.662958245626429 -6.8 -3.9999999999999996 -4.801075362242731 -5.561981460823208 -6.24456328669144 -6.814593336884906 -7.243487924448708 -7.509740483854564 -7.6 -3.9999999999999996 -4.764051981734492 -5.504888630632941 -6.200000000000001 -6.828265482620774 -7.370595549723503 -7.810511776651531 -8.134647531457997 -8.333154113253716 -8.4], ] @testset "$(shape_names[i])" for i in eachindex(shape_names) From d4ffaaab4d328a23a27cad25e3f6ac3ac78d2cae Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 4 Jan 2024 17:08:34 +0100 Subject: [PATCH 08/12] Suppress TimerOutputs warning (#329) * Suppress TimerOutputs warning * Add comments --- examples/n_body/n_body_benchmark_trixi.jl | 17 ++++++++++++++--- examples/n_body/n_body_solar_system.jl | 17 ++++++++++++++--- 2 files changed, 28 insertions(+), 6 deletions(-) diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index cb9a8d05b..6dbdf5927 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -101,7 +101,14 @@ function symplectic_euler!(velocities, coordinates, semi) end # One RHS evaluation is so fast that timers make it multiple times slower. -TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) +# Disabling timers throws a warning, which we suppress here in order to make the tests pass. +# For some reason, this only works with a file and not with devnull. See issue #332. +filename = tempname() +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + end +end @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi)) @@ -112,5 +119,9 @@ end @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi)) -# Enable timers again. -TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) +# Enable timers again +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + end +end diff --git a/examples/n_body/n_body_solar_system.jl b/examples/n_body/n_body_solar_system.jl index f063d0599..a14114fb0 100644 --- a/examples/n_body/n_body_solar_system.jl +++ b/examples/n_body/n_body_solar_system.jl @@ -44,7 +44,14 @@ saving_callback = SolutionSavingCallback(dt=10day) callbacks = CallbackSet(info_callback, saving_callback) # One RHS evaluation is so fast that timers make it multiple times slower. -TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) +# Disabling timers throws a warning, which we suppress here in order to make the tests pass. +# For some reason, this only works with a file and not with devnull. See issue #332. +filename = tempname() +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + end +end sol = solve(ode, SymplecticEuler(), dt=1.0e5, @@ -53,5 +60,9 @@ sol = solve(ode, SymplecticEuler(), @printf("%.9e\n", energy(ode.u0.x..., particle_system, semi)) @printf("%.9e\n", energy(sol.u[end].x..., particle_system, semi)) -# Enable timers again. -TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) +# Enable timers again +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + end +end From 7eb87310e686ea8fe25f1cca4e2b4f2c3a96dbaf Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 4 Jan 2024 17:08:22 +0000 Subject: [PATCH 09/12] Bump crate-ci/typos from 1.16.23 to 1.16.26 (#328) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.16.23 to 1.16.26. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.16.23...v1.16.26) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 366cb1183..a780e9751 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.16.23 + uses: crate-ci/typos@v1.16.26 From 88d287d00920ce6b88e3764e1850eefa07e5a566 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 4 Jan 2024 18:47:39 +0100 Subject: [PATCH 10/12] Add function barrier to make Adami pressure extrapolation work on macOS ARM (#330) --- .../dummy_particles/dummy_particles.jl | 34 ++++++++++++------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index a99387cdf..73432d9c7 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -364,21 +364,31 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, end @trixi_timeit timer() "inverse state equation" @threaded for particle in eachparticle(system) - # The summation is only over fluid particles, thus the volume stays zero when a boundary - # particle isn't surrounded by fluid particles. - # Check the volume to avoid NaNs in pressure and velocity. - if volume[particle] > eps() - pressure[particle] /= volume[particle] - - # To impose no-slip condition - compute_wall_velocity!(viscosity, system, system_coords, particle) - end - - # Apply inverse state equation to compute density (not used with EDAC) - inverse_state_equation!(density, state_equation, pressure, particle) + compute_adami_density!(boundary_model, system, system_coords, particle) end end +# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). +# Otherwise, `@threaded` does not work here with Julia ARM on macOS. +# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. +function compute_adami_density!(boundary_model, system, system_coords, particle) + (; pressure, state_equation, cache, viscosity) = boundary_model + (; volume, density) = cache + + # The summation is only over fluid particles, thus the volume stays zero when a boundary + # particle isn't surrounded by fluid particles. + # Check the volume to avoid NaNs in pressure and velocity. + if volume[particle] > eps() + pressure[particle] /= volume[particle] + + # To impose no-slip condition + compute_wall_velocity!(viscosity, system, system_coords, particle) + end + + # Apply inverse state equation to compute density (not used with EDAC) + inverse_state_equation!(density, state_equation, pressure, particle) +end + function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZeroing}, system, v, u, v_ode, u_ode, semi) # No pressure update needed with `PressureMirroring` and `PressureZeroing`. From f789ffb62c7acc3ecaf75e3d204329addb208b44 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 5 Jan 2024 11:34:57 +0100 Subject: [PATCH 11/12] Mixed kernel gradient (#261) * add gradient correction stuff * impl * up * finish impl * update * fix testcases * fix * fix tests * format * working * format * format * fix * working trixi_include * fix * cleanup * format * fix * fix tests * fix * ups * typo * up * fix * format * update fix errors * fix * update * up * fix tests * fix tests * np * remove comment * revert some changes * up * up * update * format * cleanup * revert * revert * implement correction on wall * fix tests * fix pressure_acceleration for summation * format * revert * fix * fixes * fix * fix pressure_acceleration dispatches * typo * fix * rename * some review fixes * review * more fixes * add solid system type * review * fix tlsph call * fix comment * review * review * format * review * move * fix doc * review fix * more review * review fix * format * change parameter list * reorder params * format * fix * fix * fix * remove dispatched functions * format * fix comments * fix * remove system * format * change parameter list * review * unused variable * break * edit todo * remove breaks * format * rename grad_kernel to W_a --- examples/fluid/dam_break_2d.jl | 13 +- examples/fluid/dam_break_2d_corrections.jl | 76 ++++-- src/TrixiParticles.jl | 5 +- src/general/corrections.jl | 251 +++++++++++++++--- src/general/density_calculators.jl | 15 -- src/general/general.jl | 10 + src/general/smoothing_kernels.jl | 33 ++- src/general/system.jl | 8 +- .../dummy_particles/dummy_particles.jl | 209 ++++++++++++--- .../monaghan_kajtar/monaghan_kajtar.jl | 12 +- src/schemes/boundary/system.jl | 18 +- src/schemes/fluid/fluid.jl | 4 - .../fluid/weakly_compressible_sph/rhs.jl | 72 +++-- .../fluid/weakly_compressible_sph/system.jl | 62 ++++- src/schemes/solid/total_lagrangian_sph/rhs.jl | 16 +- .../solid/total_lagrangian_sph/system.jl | 32 +-- test/general/smoothing_kernels.jl | 4 +- .../fluid/weakly_compressible_sph/rhs.jl | 12 +- test/systems/wcsph_system.jl | 6 +- 19 files changed, 653 insertions(+), 205 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 6a2551246..4f3b7d5fa 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -13,7 +13,7 @@ using OrdinaryDiffEq fluid_particle_spacing = 0.02 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model -boundary_layers = 4 +boundary_layers = 3 spacing_ratio = 1 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio @@ -31,7 +31,8 @@ fluid_density = 1000.0 atmospheric_pressure = 100000.0 sound_speed = 20 * sqrt(gravity * initial_fluid_size[2]) state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure, - background_pressure=atmospheric_pressure) + background_pressure=atmospheric_pressure, + clip_negative_pressure=false) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, @@ -45,6 +46,7 @@ smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, @@ -58,7 +60,8 @@ boundary_density_calculator = AdamiPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + correction=nothing) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) @@ -68,7 +71,7 @@ semi = Semidiscretization(fluid_system, boundary_system, neighborhood_search=GridNeighborhoodSearch) ode = semidiscretize(semi, tspan) -info_callback = InfoCallback(interval=100) +info_callback = InfoCallback(interval=250) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") use_reinit = false @@ -87,5 +90,5 @@ callbacks = CallbackSet(info_callback, saving_callback, density_reinit_cb) sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-5, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing + dtmax=1e-3, # Limit stepsize to prevent crashing save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_2d_corrections.jl b/examples/fluid/dam_break_2d_corrections.jl index ad697d454..bd2b790af 100644 --- a/examples/fluid/dam_break_2d_corrections.jl +++ b/examples/fluid/dam_break_2d_corrections.jl @@ -3,19 +3,35 @@ using TrixiParticles fluid_density = 1000.0 particle_spacing = 0.05 -smoothing_length = 1.15 * particle_spacing -boundary_density_calculator = SummationDensity() - -gravity = 9.81 -tspan = (0.0, 5.7 / sqrt(gravity)) +tspan = (0.0, 5.7 / sqrt(9.81)) correction_dict = Dict( - "no_correction" => Nothing(), + "no_correction" => nothing, "shepard_kernel_correction" => ShepardKernelCorrection(), "akinci_free_surf_correction" => AkinciFreeSurfaceCorrection(fluid_density), - "kernel_gradient_summation_correction" => KernelGradientCorrection(), - "kernel_gradient_continuity_correction" => KernelGradientCorrection(), + "kernel_gradient_summation_correction" => KernelCorrection(), + "kernel_gradient_continuity_correction" => KernelCorrection(), + "blended_gradient_summation_correction" => BlendedGradientCorrection(0.5), + "blended_gradient_continuity_correction" => BlendedGradientCorrection(0.2), + "gradient_summation_correction" => GradientCorrection(), + "mixed_kernel_gradient_summation_correction" => MixedKernelGradientCorrection(), + # "gradient_continuity_correction" => GradientCorrection(), + # "mixed_kernel_gradient_continuity_correction" => MixedKernelGradientCorrection(), +) + +smoothing_length_dict = Dict( + "no_correction" => 3.0 * particle_spacing, + "shepard_kernel_correction" => 3.0 * particle_spacing, + "akinci_free_surf_correction" => 3.0 * particle_spacing, + "kernel_gradient_summation_correction" => 4.0 * particle_spacing, + "kernel_gradient_continuity_correction" => 3.5 * particle_spacing, + "blended_gradient_summation_correction" => 3.0 * particle_spacing, + "blended_gradient_continuity_correction" => 4.0 * particle_spacing, + "gradient_summation_correction" => 3.5 * particle_spacing, + "mixed_kernel_gradient_summation_correction" => 3.5 * particle_spacing, + "gradient_continuity_correction" => 4.5 * particle_spacing, + "mixed_kernel_gradient_continuity_correction" => 4.0 * particle_spacing, ) density_calculator_dict = Dict( @@ -24,23 +40,42 @@ density_calculator_dict = Dict( "akinci_free_surf_correction" => SummationDensity(), "kernel_gradient_summation_correction" => SummationDensity(), "kernel_gradient_continuity_correction" => ContinuityDensity(), + "blended_gradient_summation_correction" => SummationDensity(), + "blended_gradient_continuity_correction" => ContinuityDensity(), + "gradient_summation_correction" => SummationDensity(), + "gradient_continuity_correction" => ContinuityDensity(), + "mixed_kernel_gradient_summation_correction" => SummationDensity(), + "mixed_kernel_gradient_continuity_correction" => ContinuityDensity(), +) + +smoothing_kernel_dict = Dict( + "no_correction" => WendlandC2Kernel{2}(), + "shepard_kernel_correction" => WendlandC2Kernel{2}(), + "akinci_free_surf_correction" => WendlandC2Kernel{2}(), + "kernel_gradient_summation_correction" => WendlandC6Kernel{2}(), + "kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(), + "blended_gradient_summation_correction" => WendlandC2Kernel{2}(), + "blended_gradient_continuity_correction" => WendlandC6Kernel{2}(), + "gradient_summation_correction" => WendlandC6Kernel{2}(), + "gradient_continuity_correction" => WendlandC6Kernel{2}(), + "mixed_kernel_gradient_summation_correction" => WendlandC6Kernel{2}(), + "mixed_kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(), ) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), - fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length, + fluid_particle_spacing=particle_spacing, + smoothing_length=3.0 * particle_spacing, boundary_density_calculator=ContinuityDensity(), fluid_density_calculator=ContinuityDensity(), - correction=Nothing(), use_reinit=true, - prefix="continuity_reinit", tspan=tspan) - -# Clip negative pressure to be able to use `SummationDensity` -state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure, - background_pressure=atmospheric_pressure, - clip_negative_pressure=true) + correction=nothing, use_reinit=true, + prefix="continuity_reinit", tspan=tspan, + fluid_density=fluid_density, density_diffusion=nothing) for correction_name in keys(correction_dict) local fluid_density_calculator = density_calculator_dict[correction_name] local correction = correction_dict[correction_name] + local smoothing_kernel = smoothing_kernel_dict[correction_name] + local smoothing_length = smoothing_length_dict[correction_name] println("="^100) println("fluid/dam_break_2d.jl with ", correction_name) @@ -48,9 +83,12 @@ for correction_name in keys(correction_dict) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length, - boundary_density_calculator=boundary_density_calculator, + boundary_density_calculator=SummationDensity(), fluid_density_calculator=fluid_density_calculator, correction=correction, use_reinit=false, - state_equation=state_equation, - prefix="$(correction_name)", tspan=tspan) + clip_negative_pressure=(fluid_density_calculator isa SummationDensity), + smoothing_kernel=smoothing_kernel, + prefix="$(correction_name)", tspan=tspan, + fluid_density=fluid_density, density_diffusion=Nothing(), + boundary_layers=5) end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c9bff5557..87789eeb1 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -5,7 +5,7 @@ using Reexport: @reexport using Dates using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect using FastPow: @fastpow -using LinearAlgebra: norm, dot, I, tr +using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using Morton: cartesian2morton using MuladdMacro: @muladd using Polyester: Polyester, @batch @@ -55,7 +55,8 @@ export examples_dir, trixi_include export trixi2vtk export RectangularTank, RectangularShape, SphereShape export VoxelSphere, RoundSphere, reset_wall! -export ShepardKernelCorrection, KernelGradientCorrection, AkinciFreeSurfaceCorrection +export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, + GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection export nparticles end # module diff --git a/src/general/corrections.jl b/src/general/corrections.jl index 595a61f0e..21fc8f89d 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -1,5 +1,4 @@ # Sorted in order of computational cost - @doc raw""" AkinciFreeSurfaceCorrection(rho0) @@ -11,6 +10,16 @@ near free surfaces to counter this effect. It's important to note that this correlation is unphysical and serves as an approximation. The computation time added by this method is about 2--3%. +Mathematically the idea is quite simple. If we have an SPH particle in the middle of a volume +at rest, its density will be identical to the rest density ``\rho_0``. If we now consider an SPH +particle at a free surface at rest, it will have neighbors missing in the direction normal to +the surface, which will result in a lower density. If we calculate the correction factor +```math +k = \rho_0/\rho_\text{mean}, +``` +this value will be about ~1.5 for particles at the free surface and can then be used to increase +the pressure and viscosity accordingly. + # Arguments - `rho0`: Rest density. @@ -71,19 +80,19 @@ to an improvement, especially at free surfaces. [doi: 10.1016/S0045-7825(99)00051-1](https://doi.org/10.1016/S0045-7825(99)00051-1) - Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". - In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127-1148. + In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127--1148. [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927) - Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". - In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159-193. + In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159--193. [doi:10.1016/S0045-7825(96)01082-1](https://doi.org/10.1016/S0045-7825(96)01082-1) """ struct ShepardKernelCorrection end @doc raw""" - KernelGradientCorrection() + KernelCorrection() -Kernel gradient correction uses Shepard interpolation to obtain a 0-th order accurate result, which +Kernel correction uses Shepard interpolation to obtain a 0-th order accurate result, which was first proposed by Li et al. This can be further extended to obtain a kernel corrected gradient as shown by Basa et al. @@ -93,7 +102,7 @@ c(x) = \sum_{b=1} V_b W_b(x) ``` The gradient of corrected kernel is determined by ```math -\nabla \tilde{W}_{b}(r) =\frac{\nabla W_{b}(r) - \gamma(r)}{\sum_{b=1} V_b W_b(r)} , \quad \text{where} \quad +\nabla \tilde{W}_{b}(r) =\frac{\nabla W_{b}(r) - W_b(r) \gamma(r)}{\sum_{b=1} V_b W_b(r)} , \quad \text{where} \quad \gamma(r) = \frac{\sum_{b=1} V_b \nabla W_b(r)}{\sum_{b=1} V_b W_b(r)}. ``` @@ -113,30 +122,64 @@ This correction can be applied with [`SummationDensity`](@ref) and [doi: 10.1016/S0045-7825(99)00051-1](https://doi.org/10.1016/S0045-7825(99)00051-1) - Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". - In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127-1148. + In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127--1148. [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927) - Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159-193. [doi:10.1016/S0045-7825(96)01082-1](https://doi.org/10.1016/S0045-7825(96)01082-1) """ -struct KernelGradientCorrection end +struct KernelCorrection end -function kernel_correction_coefficient(system, particle) +@doc raw""" + MixedKernelGradientCorrection() + +Combines [`GradientCorrection`](@ref) and [`KernelCorrection`](@ref), +which results in a 1st-order-accurate SPH method. + +# Notes: +- Stability issues, especially when particles separate into small clusters. +- Doubles the computational effort. + +## References: +- J. Bonet, T.-S.L. Lok. + "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". + In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97--115. + [doi: 10.1016/S0045-7825(99)00051-1](https://doi.org/10.1016/S0045-7825(99)00051-1) +- Mihai Basa, Nathan Quinlan, Martin Lastiwka. + "Robustness and accuracy of SPH formulations for viscous flow". + In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127--1148. + [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927) +""" +struct MixedKernelGradientCorrection end + +function kernel_correction_coefficient(system::FluidSystem, particle) return system.cache.kernel_correction_coefficient[particle] end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - density_calculator, correction) +function kernel_correction_coefficient(system::BoundarySystem, particle) + return system.boundary_model.cache.kernel_correction_coefficient[particle] +end + +function compute_correction_values!(system, correction, v, u, v_ode, u_ode, semi, + density_calculator) return system end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - ::SummationDensity, ::ShepardKernelCorrection) +function compute_correction_values!(system, ::ShepardKernelCorrection, v, u, v_ode, u_ode, + semi, + ::SummationDensity) return compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, system.cache.kernel_correction_coefficient) end +function compute_correction_values!(system::BoundarySystem, ::ShepardKernelCorrection, v, u, + v_ode, u_ode, semi, + ::SummationDensity) + return compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, + system.boundary_model.cache.kernel_correction_coefficient) +end + function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, kernel_correction_coefficient) set_zero!(kernel_correction_coefficient) @@ -153,8 +196,9 @@ function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, # Loop over all pairs of particles and neighbors within the kernel cutoff for_particle_neighbor(system, neighbor_system, system_coords, - neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance + neighbor_coords, neighborhood_search, + particles=eachparticle(system)) do particle, neighbor, + pos_diff, distance rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) m_b = hydrodynamic_mass(neighbor_system, neighbor) volume = m_b / rho_b @@ -167,16 +211,43 @@ function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, return kernel_correction_coefficient end -function dw_gamma(system, particle) +function dw_gamma(system::FluidSystem, particle) return extract_svector(system.cache.dw_gamma, system, particle) end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - ::Union{SummationDensity, ContinuityDensity}, - ::KernelGradientCorrection) - (; cache) = system - (; kernel_correction_coefficient, dw_gamma) = cache +function dw_gamma(system::BoundarySystem, particle) + return extract_svector(system.boundary_model.cache.dw_gamma, system, particle) +end + +function compute_correction_values!(system::FluidSystem, + correction::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, + v_ode, u_ode, semi, + density_calculator) + compute_correction_values!(system, correction, + v, u, v_ode, u_ode, semi, + density_calculator, + system.cache.kernel_correction_coefficient, + system.cache.dw_gamma) +end +function compute_correction_values!(system::BoundarySystem, + correction::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, + v_ode, u_ode, semi, + density_calculator) + compute_correction_values!(system, correction, v, u, v_ode, u_ode, semi, + density_calculator, + system.boundary_model.cache.kernel_correction_coefficient, + system.boundary_model.cache.dw_gamma) +end + +function compute_correction_values!(system, + ::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, v_ode, + u_ode, semi, + density_calculator, kernel_correction_coefficient, + dw_gamma) set_zero!(kernel_correction_coefficient) set_zero!(dw_gamma) @@ -193,7 +264,9 @@ function compute_correction_values!(system, v, u, v_ode, u_ode, semi, # Loop over all pairs of particles and neighbors within the kernel cutoff for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance + neighborhood_search, + particles=eachparticle(system)) do particle, neighbor, + pos_diff, distance rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) m_b = hydrodynamic_mass(neighbor_system, neighbor) volume = m_b / rho_b @@ -250,21 +323,46 @@ the correction matrix $\bm{L}_a$ is evaluated explicitly as !!! note - Stability issues arise, especially when particles separate into small clusters. - Doubles the computational effort. +- Better stability with smoother smoothing Kernels with larger support, e.g. [`SchoenbergQuinticSplineKernel`](@ref) or [`WendlandC6Kernel`](@ref). +- Set `dt_max =< 1e-3` for stability. ## References: - J. Bonet, T.-S.L. Lok. "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". - In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97-115. + In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97--115. [doi: 10.1016/S0045-7825(99)00051-1](https://doi.org/10.1016/S0045-7825(99)00051-1) - Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". - In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127-1148. + In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127--1148. [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927) """ struct GradientCorrection end -function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, system, - coordinates, density_fun) +@doc raw""" + BlendedGradientCorrection() + +Calculate a blended gradient to reduce the stability issues of the [`GradientCorrection`](@ref). + +This calculates the following, +```math +\tilde\nabla A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i +``` +with ``0 \leq \lambda \leq 1`` being the blending factor. + +# Arguments +- `blending_factor`: Blending factor between corrected and regular SPH gradient. +""" +struct BlendedGradientCorrection{ELTYPE <: Real} + blending_factor::ELTYPE + + function BlendedGradientCorrection(blending_factor) + return new{eltype(blending_factor)}(blending_factor) + end +end + +# Called only by DensityDiffusion and TLSPH +function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, + system, coordinates, density_fun) (; mass) = system set_zero!(corr_matrix) @@ -275,12 +373,12 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s neighborhood_search; particles=eachparticle(system)) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0. - distance < sqrt(eps()) && return - volume = mass[neighbor] / density_fun(neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + + iszero(grad_kernel) && return + result = volume * grad_kernel * pos_diff' @inbounds for j in 1:ndims(system), i in 1:ndims(system) @@ -288,20 +386,95 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s end end - @threaded for particle in eachparticle(system) - L = extract_smatrix(corr_matrix, system, particle) - result = inv(L) + correction_matrix_inversion_step!(corr_matrix, system) - if any(isinf.(result)) || any(isnan.(result)) - # TODO How do we handle singularities correctly? - # See https://github.com/trixi-framework/TrixiParticles.jl/issues/273 - result = one(L) - end + return corr_matrix +end - @inbounds for j in 1:ndims(system), i in 1:ndims(system) - corr_matrix[i, j, particle] = result[i, j] +function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, + coordinates, v_ode, u_ode, semi, + correction, smoothing_length, smoothing_kernel) + set_zero!(corr_matrix) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + @trixi_timeit timer() "compute correction matrix" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + neighborhood_search = neighborhood_searches(system, neighbor_system, semi) + + for_particle_neighbor(system, neighbor_system, coordinates, neighbor_coords, + neighborhood_search; + particles=eachparticle(system)) do particle, + neighbor, + pos_diff, + distance + volume = hydrodynamic_mass(neighbor_system, neighbor) / + particle_density(v_neighbor_system, neighbor_system, neighbor) + + function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, + smoothing_length, system, particle) + return smoothing_kernel_grad(system, pos_diff, distance) + end + + # Compute gradient of corrected kernel + function compute_grad_kernel(correction::MixedKernelGradientCorrection, + smoothing_kernel, pos_diff, distance, + smoothing_length, system, particle) + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length, KernelCorrection(), system, + particle) + end + + grad_kernel = compute_grad_kernel(correction, smoothing_kernel, pos_diff, + distance, smoothing_length, system, particle) + + iszero(grad_kernel) && return + + L = volume * grad_kernel * pos_diff' + + # pos_diff is always x_a - x_b hence * -1 to switch the order to x_b - x_a + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] -= L[i, j] + end end end + correction_matrix_inversion_step!(corr_matrix, system) + return corr_matrix end + +function correction_matrix_inversion_step!(corr_matrix, system) + @threaded for particle in eachparticle(system) + L = extract_smatrix(corr_matrix, system, particle) + norm_ = norm(L) + + # The norm value is quasi-zero, so there are probably no neighbors for this particle + if norm_ < sqrt(eps()) + # The correction matrix is set to an identity matrix, which effectively disables + # the correction for this particle. + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] = 0.0 + end + @inbounds for i in 1:ndims(system) + corr_matrix[i, i, particle] = 1.0 + end + continue + end + + det_ = abs(det(L)) + @fastmath if det_ < 1e-6 * norm_ + L_inv = pinv(L) + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] = L_inv[i, j] + end + else + L_inv = inv(L) + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] = L_inv[i, j] + end + end + end +end diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 4383d55e3..732a5eb2a 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -35,21 +35,6 @@ end return v[end, particle] end -# *Note* that these functions are intended to internally set the density for buffer particles -# and density correction. It cannot be used to set up an initial condition, -# as the particle density depends on the particle positions. -@inline function set_particle_density(particle, v, system, density) - set_particle_density(particle, v, system.density_calculator, system, density) -end - -@inline function set_particle_density(particle, v, ::SummationDensity, system, density) - system.cache.density[particle] = density -end - -@inline function set_particle_density(particle, v, ::ContinuityDensity, system, density) - v[end, particle] = density -end - function summation_density!(system, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) diff --git a/src/general/general.jl b/src/general/general.jl index e5cf8bdd4..a5c3f8fc5 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -1,3 +1,13 @@ +abstract type System{NDIMS} end + +abstract type FluidSystem{NDIMS} <: System{NDIMS} end + +abstract type SolidSystem{NDIMS} <: System{NDIMS} end + +abstract type BoundarySystem{NDIMS} <: System{NDIMS} end + +timer_name(::FluidSystem) = "fluid" + @inline function set_zero!(du) du .= zero(eltype(du)) diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index c71869ee5..c2f48ebba 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -52,6 +52,8 @@ abstract type SmoothingKernel{NDIMS} end @inline Base.ndims(::SmoothingKernel{NDIMS}) where {NDIMS} = NDIMS @inline function kernel_grad(kernel, pos_diff, distance, h) + distance < sqrt(eps()) && return zero(pos_diff) + return kernel_deriv(kernel, distance, h) / distance * pos_diff end @@ -60,12 +62,37 @@ end return kernel_grad(kernel, pos_diff, distance, h) end -@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, - ::KernelGradientCorrection, system, particle) - return (kernel_grad(kernel, pos_diff, distance, h) .- dw_gamma(system, particle)) / +@inline function corrected_kernel_grad(kernel_, pos_diff, distance, h, ::KernelCorrection, + system, particle) + return (kernel_grad(kernel_, pos_diff, distance, h) .- + kernel(kernel_, distance, h) * dw_gamma(system, particle)) / kernel_correction_coefficient(system, particle) end +@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, + corr::BlendedGradientCorrection, system, + particle) + (; blending_factor) = corr + + grad = kernel_grad(kernel, pos_diff, distance, h) + return (1 - blending_factor) * grad + + blending_factor * correction_matrix(system, particle) * grad +end + +@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, + ::GradientCorrection, system, particle) + grad = kernel_grad(kernel, pos_diff, distance, h) + return correction_matrix(system, particle) * grad +end + +@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, + ::MixedKernelGradientCorrection, system, + particle) + grad = corrected_kernel_grad(kernel, pos_diff, distance, h, KernelCorrection(), + system, particle) + return correction_matrix(system, particle) * grad +end + @doc raw""" GaussianKernel{NDIMS}() diff --git a/src/general/system.jl b/src/general/system.jl index a65796a66..b27b229e3 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -1,5 +1,3 @@ -abstract type System{NDIMS} end - initialize!(system, neighborhood_search) = system @inline Base.ndims(::System{NDIMS}) where {NDIMS} = NDIMS @@ -80,6 +78,12 @@ end return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) end +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance) + (; smoothing_kernel, smoothing_length) = system.boundary_model + + return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) +end + # This is dispatched for some system types @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 73432d9c7..4ca85c33e 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -74,7 +74,7 @@ We provide five options to compute the boundary density and pressure, determined In: Computers, Materials and Continua 5 (2007), pages 173-184. [doi: 10.3970/cmc.2007.005.173](https://doi.org/10.3970/cmc.2007.005.173) """ -struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, C} +struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, COR, C} pressure :: Vector{ELTYPE} hydrodynamic_mass :: Vector{ELTYPE} state_equation :: SE @@ -82,37 +82,34 @@ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, C} smoothing_kernel :: K smoothing_length :: ELTYPE viscosity :: V + correction :: COR cache :: C function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=NoViscosity(), - state_equation=nothing) + state_equation=nothing, correction=nothing) pressure = initial_boundary_pressure(initial_density, density_calculator, state_equation) + NDIMS = ndims(smoothing_kernel) n_particles = length(initial_density) - cache = (; create_cache_model(viscosity, n_particles, ndims(smoothing_kernel))..., + cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., create_cache_model(initial_density, density_calculator)...) + cache = (; + create_cache_model(correction, initial_density, NDIMS, + n_particles)..., cache...) new{typeof(density_calculator), eltype(initial_density), typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), - typeof(cache)}(pressure, hydrodynamic_mass, state_equation, density_calculator, - smoothing_kernel, smoothing_length, viscosity, cache) + typeof(correction), typeof(cache)}(pressure, hydrodynamic_mass, state_equation, + density_calculator, + smoothing_kernel, smoothing_length, + viscosity, correction, cache) end end -function Base.show(io::IO, model::BoundaryModelDummyParticles) - @nospecialize model # reduce precompilation time - - print(io, "BoundaryModelDummyParticles(") - print(io, model.density_calculator |> typeof |> nameof) - print(io, ", ") - print(io, model.viscosity |> typeof |> nameof) - print(io, ")") -end - @doc raw""" AdamiPressureExtrapolation() @@ -188,39 +185,27 @@ reference pressure (the corresponding pressure to the reference density by the s """ struct PressureZeroing end -# For most density calculators, the pressure is updated in every step -initial_boundary_pressure(initial_density, density_calculator, _) = similar(initial_density) -# Pressure mirroring does not use the pressure, so we set it to zero for the visualization -initial_boundary_pressure(initial_density, ::PressureMirroring, _) = zero(initial_density) +create_cache_model(correction, density, NDIMS, nparticles) = (;) -# For pressure zeroing, set the pressure to the reference pressure (zero with free surfaces) -function initial_boundary_pressure(initial_density, ::PressureZeroing, state_equation) - return state_equation.(initial_density) +function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density)) end -# With EDAC, just use zero pressure -function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing) - return zero(initial_density) +function create_cache_model(::KernelCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density), dw_gamma) end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, smoothing_length, - grad_kernel, - boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, - fluid_density_calculator) - - # Use `p_a` as pressure for both particles with `PressureMirroring` - return pressure_acceleration(pressure_correction, m_b, p_a, p_a, rho_a, rho_b, - grad_kernel, fluid_density_calculator) +function create_cache_model(::Union{GradientCorrection, BlendedGradientCorrection}, density, + NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; correction_matrix) end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, smoothing_length, - grad_kernel, - boundary_model::BoundaryModelDummyParticles, - fluid_density_calculator) - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, fluid_density_calculator) +function create_cache_model(::MixedKernelGradientCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end function create_cache_model(initial_density, @@ -264,6 +249,90 @@ function reset_cache!(cache, viscosity::ViscosityAdami) return cache end +function Base.show(io::IO, model::BoundaryModelDummyParticles) + @nospecialize model # reduce precompilation time + + print(io, "BoundaryModelDummyParticles(") + print(io, model.density_calculator |> typeof |> nameof) + print(io, ", ") + print(io, model.viscosity |> typeof |> nameof) + print(io, ")") +end + +# For most density calculators, the pressure is updated in every step +initial_boundary_pressure(initial_density, density_calculator, _) = similar(initial_density) +# Pressure mirroring does not use the pressure, so we set it to zero for the visualization +initial_boundary_pressure(initial_density, ::PressureMirroring, _) = zero(initial_density) + +# For pressure zeroing, set the pressure to the reference pressure (zero with free surfaces) +function initial_boundary_pressure(initial_density, ::PressureZeroing, state_equation) + return state_equation.(initial_density) +end + +# With EDAC, just use zero pressure +function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing) + return zero(initial_density) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, + fluid_density_calculator, + correction) + + # Use `p_a` as pressure for both particles with `PressureMirroring` + return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_a, rho_a, rho_b, + grad_kernel, fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles, + fluid_density_calculator, + correction) + return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, + grad_kernel, fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, W_a, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, + fluid_density_calculator, + correction::Union{KernelCorrection, + GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}) + W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor) + + # Use `p_a` as pressure for both particles with `PressureMirroring` + return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_a, + rho_a, rho_b, W_a, W_b, + fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, W_a, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles, + fluid_density_calculator, + correction::Union{KernelCorrection, + GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}) + W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor) + + return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, W_a, W_b, + fluid_density_calculator) +end + @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle) return particle_density(v, model.density_calculator, model, particle) end @@ -302,15 +371,55 @@ end @inline function update_pressure!(boundary_model::BoundaryModelDummyParticles, system, v, u, v_ode, u_ode, semi) - (; density_calculator) = boundary_model + (; density_calculator, correction) = boundary_model + + compute_correction_values!(system, + correction, v, u, v_ode, u_ode, semi, density_calculator) + + compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode, + semi) + + # `kernel_correct_density!` only performed for `SummationDensity` + kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction, + density_calculator) compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) return boundary_model end -function compute_density!(boundary_model, ::SummationDensity, - system, v, u, v_ode, u_ode, semi) +function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, + correction, density_calculator) + return boundary_model +end + +function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, + corr::ShepardKernelCorrection, ::SummationDensity) + boundary_model.cache.density ./= boundary_model.cache.kernel_correction_coefficient +end + +function compute_gradient_correction_matrix!(correction, boundary_model, system, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + boundary_model, + system, u, v_ode, u_ode, semi) + (; cache, correction, smoothing_kernel, smoothing_length) = boundary_model + (; correction_matrix) = cache + + system_coords = current_coordinates(u, system) + + compute_gradient_correction_matrix!(correction_matrix, system, system_coords, + v_ode, u_ode, semi, correction, smoothing_length, + smoothing_kernel) +end + +function compute_density!(boundary_model, ::SummationDensity, system, v, u, v_ode, u_ode, + semi) (; cache) = boundary_model (; density) = cache # Density is in the cache for SummationDensity @@ -487,3 +596,15 @@ end # The density is constant when using EDAC return density end + +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, + distance, particle) + (; smoothing_kernel, smoothing_length, correction) = system.boundary_model + + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length, correction, system, particle) +end + +@inline function correction_matrix(system::BoundarySystem, particle) + extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) +end diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index e608128a6..fa476e530 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -76,11 +76,13 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff::SVector{NDIMS}, - smoothing_length, grad_kernel, - boundary_model::BoundaryModelMonaghanKajtar, - density_calculator) where {NDIMS} +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff::SVector{NDIMS}, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelMonaghanKajtar, + fluid_density_calculator, + correction) where {NDIMS} (; K, beta, boundary_particle_spacing) = boundary_model distance = norm(pos_diff) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index e3cc8c29d..a5be4ab93 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -7,7 +7,7 @@ The interaction between fluid and boundary particles is specified by the boundar For moving boundaries, a [`BoundaryMovement`](@ref) can be passed with the keyword argument `movement`. """ -struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: System{NDIMS} +struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIMS} coordinates :: Array{ELTYPE, 2} boundary_model :: BM movement :: M @@ -266,13 +266,17 @@ function viscosity_model(system::BoundarySPHSystem) end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system::BoundarySPHSystem, - density_calculator) + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, + neighbor_system::BoundarySPHSystem, + density_calculator, correction) (; boundary_model) = neighbor_system (; smoothing_length) = particle_system - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - smoothing_length, grad_kernel, boundary_model, - density_calculator) + return pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model, + density_calculator, correction) end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 844e0a91d..712a6c7dc 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -1,7 +1,3 @@ -abstract type FluidSystem{NDIMS} <: System{NDIMS} end - -timer_name(::FluidSystem) = "fluid" - @inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] function write_u0!(u0, system::FluidSystem) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 11a9e8d4e..f0461f477 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -21,9 +21,6 @@ function interact!(dv, v_particle_system, u_particle_system, for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return - rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = 0.5 * (rho_a + rho_b) @@ -42,9 +39,9 @@ function interact!(dv, v_particle_system, u_particle_system, p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) dv_pressure = pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system, - density_calculator) + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, neighbor_system, + density_calculator, correction) dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, @@ -58,6 +55,7 @@ function interact!(dv, v_particle_system, u_particle_system, # debug_array[i, particle] += dv_pressure[i] end + # TODO If variable smoothing_length is used, this should use the neighbor smoothing length continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, particle_system, neighbor_system, grad_kernel) @@ -73,13 +71,33 @@ function interact!(dv, v_particle_system, u_particle_system, end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system, - density_calculator) + rho_a, rho_b, pos_diff, distance, + grad_kernel, particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + density_calculator, + correction) + + # Without correction, the kernel gradient is symmetric, so call the symmetric + # pressure acceleration formulation corresponding to the density calculator. + return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, + grad_kernel, density_calculator) +end - # By default, just call the pressure acceleration formulation corresponding to the density calculator - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, density_calculator) +@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + W_a, particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + density_calculator, + correction::Union{KernelCorrection, + GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}) + W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor) + + # With correction, the kernel gradient is not necessarily symmetric, so call the + # asymmetric pressure acceleration formulation corresponding to the density calculator. + return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b, rho_a, + rho_b, W_a, W_b, density_calculator) end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic @@ -87,8 +105,10 @@ end # used with `ContinuityDensity` with the formulation `\rho_a * \sum m_b / \rho_b ...`. # This can also be seen in the tests for total energy conservation, which fail with the # other `pressure_acceleration` form. -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, ::ContinuityDensity) +# We assume symmetry of the kernel gradient in this formulation. See below for the +# asymmetric version. +@inline function pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, + rho_b, grad_kernel, ::ContinuityDensity) return (-m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel) * pressure_correction end @@ -97,11 +117,27 @@ end # used with `SummationDensity`. # This can also be seen in the tests for total energy conservation, which fail with the # other `pressure_acceleration` form. -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, ::SummationDensity) +# We assume symmetry of the kernel gradient in this formulation. See below for the +# asymmetric version. +@inline function pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, + rho_b, grad_kernel, ::SummationDensity) return (-m_b * (p_a / rho_a^2 + p_b / rho_b^2) * grad_kernel) * pressure_correction end +# Same as above, but not assuming symmetry of the kernel gradient. To be used with +# corrections that do not produce a symmetric kernel gradient. +@inline function pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b, rho_a, + rho_b, W_a, W_b, ::ContinuityDensity) + return -m_b / (rho_a * rho_b) * (p_a * W_a - p_b * W_b) * pressure_correction +end + +# Same as above, but not assuming symmetry of the kernel gradient. To be used with +# corrections that do not produce a symmetric kernel gradient. +@inline function pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b, rho_a, + rho_b, W_a, W_b, ::SummationDensity) + return (-m_b * (p_a / rho_a^2 * W_a - p_b / rho_b^2 * W_b)) * pressure_correction +end + # With 'SummationDensity', density is calculated in wcsph/system.jl:compute_density! @inline function continuity_equation!(dv, density_calculator::SummationDensity, v_particle_system, v_neighbor_system, @@ -111,6 +147,7 @@ end return dv end +# This formulation was chosen to be consistent with the used pressure_acceleration formulations. @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -136,6 +173,9 @@ end particle_system::WeaklyCompressibleSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem, grad_kernel) + # Density diffusion terms are all zero for distance zero + distance < sqrt(eps()) && return + (; delta) = density_diffusion (; smoothing_length, state_equation) = particle_system (; sound_speed) = state_equation diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index e65d1465a..4108466a6 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -87,12 +87,25 @@ end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) function create_cache_wcsph(::ShepardKernelCorrection, density, NDIMS, n_particles) - (; kernel_correction_coefficient=similar(density)) + return (; kernel_correction_coefficient=similar(density)) end -function create_cache_wcsph(::KernelGradientCorrection, density, NDIMS, n_particles) +function create_cache_wcsph(::KernelCorrection, density, NDIMS, n_particles) dw_gamma = Array{Float64}(undef, NDIMS, n_particles) - (; kernel_correction_coefficient=similar(density), dw_gamma) + return (; kernel_correction_coefficient=similar(density), dw_gamma) +end + +function create_cache_wcsph(::Union{GradientCorrection, BlendedGradientCorrection}, density, + NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; correction_matrix) +end + +function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + + return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end function create_cache_wcsph(n_particles, ELTYPE, ::SummationDensity) @@ -162,7 +175,7 @@ initialize!(system::WeaklyCompressibleSPHSystem, neighborhood_search) = system function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, density_diffusion) = system + (; density_calculator, density_diffusion, correction) = system compute_density!(system, u, u_ode, semi, density_calculator) @@ -188,8 +201,11 @@ end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, correction) = system - compute_correction_values!(system, v, u, v_ode, u_ode, semi, - density_calculator, correction) + compute_correction_values!(system, + correction, v, u, v_ode, u_ode, semi, density_calculator) + + compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) + # `kernel_correct_density!` only performed for `SummationDensity` kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) @@ -198,17 +214,37 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end -function kernel_correct_density!(system, v, u, v_ode, u_ode, semi, - correction, density_calculator) +function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, + semi, correction, density_calculator) return system end -function kernel_correct_density!(system, v, u, v_ode, u_ode, semi, - ::Union{ShepardKernelCorrection, KernelGradientCorrection}, - ::SummationDensity) +function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, + semi, corr::ShepardKernelCorrection, ::SummationDensity) system.cache.density ./= system.cache.kernel_correction_coefficient end +function compute_gradient_correction_matrix!(correction, + system::WeaklyCompressibleSPHSystem, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + system::WeaklyCompressibleSPHSystem, u, + v_ode, u_ode, semi) + (; cache, correction, smoothing_kernel, smoothing_length) = system + (; correction_matrix) = cache + + system_coords = current_coordinates(u, system) + + compute_gradient_correction_matrix!(correction_matrix, system, system_coords, + v_ode, u_ode, semi, correction, smoothing_length, + smoothing_kernel) +end + function reinit_density!(vu_ode, semi) v_ode, u_ode = vu_ode.x @@ -313,3 +349,7 @@ end system.smoothing_length, system.correction, system, particle) end + +@inline function correction_matrix(system::WeaklyCompressibleSPHSystem, particle) + extract_smatrix(system.cache.correction_matrix, system, particle) +end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index bd9505505..f5c17dd68 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -61,7 +61,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem) (; boundary_model) = particle_system - (; density_calculator, state_equation, viscosity) = neighbor_system + (; density_calculator, state_equation, viscosity, correction) = neighbor_system (; sound_speed) = state_equation system_coords = current_coordinates(u_particle_system, particle_system) @@ -91,10 +91,9 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) # use `m_a` to get the same viscosity as for the fluid-solid direction. - dv_viscosity = viscosity(neighbor_system, particle_system, - v_neighbor_system, v_particle_system, - neighbor, particle, - pos_diff, distance, sound_speed, m_b, m_a, rho_mean) + dv_viscosity = viscosity(neighbor_system, particle_system, v_neighbor_system, + v_particle_system, neighbor, particle, pos_diff, distance, + sound_speed, m_b, m_a, rho_mean) # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles # corresponding to the chosen boundary model. @@ -104,9 +103,10 @@ function interact!(dv, v_particle_system, u_particle_system, # Boundary forces # Note: neighbor and particle pressure are switched in this call # and `pressure_correction` is set to `1.0` (no correction) - dv_boundary = pressure_acceleration(1.0, m_b, p_b, p_a, rho_b, rho_a, pos_diff, - neighbor_system.smoothing_length, - grad_kernel, boundary_model, density_calculator) + dv_boundary = pressure_acceleration(1.0, m_a, p_b, p_a, + rho_b, rho_a, pos_diff, distance, grad_kernel, + neighbor_system, neighbor, particle_system, + density_calculator, correction) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 013048332..14a4b37ab 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -75,7 +75,7 @@ The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGa In: International Journal for Numerical Methods in Engineering 48 (2000), pages 1359–1400. [doi: 10.1002/1097-0207](https://doi.org/10.1002/1097-0207) """ -struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIMS} +struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} initial_coordinates :: Array{ELTYPE, 2} # [dimension, particle] current_coordinates :: Array{ELTYPE, 2} # [dimension, particle] @@ -94,7 +94,6 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM acceleration :: SVector{NDIMS, ELTYPE} boundary_model :: BM penalty_force :: PF - function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio, boundary_model; @@ -130,17 +129,15 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) lame_mu = 0.5 * young_modulus / (1 + poisson_ratio) - return new{typeof(boundary_model), - NDIMS, ELTYPE, + return new{typeof(boundary_model), NDIMS, ELTYPE, typeof(smoothing_kernel), typeof(penalty_force)}(initial_condition, initial_coordinates, - current_coordinates, mass, - correction_matrix, pk1_corrected, - deformation_grad, material_density, + current_coordinates, mass, correction_matrix, + pk1_corrected, deformation_grad, material_density, n_moving_particles, young_modulus, poisson_ratio, - lame_lambda, lame_mu, - smoothing_kernel, smoothing_length, - acceleration_, boundary_model, penalty_force) + lame_lambda, lame_mu, smoothing_kernel, + smoothing_length, acceleration_, boundary_model, + penalty_force) end end @@ -234,6 +231,7 @@ end @inline function correction_matrix(system, particle) extract_smatrix(system.correction_matrix, system, particle) end + @inline function deformation_gradient(system, particle) extract_smatrix(system.deformation_grad, system, particle) end @@ -412,16 +410,18 @@ function viscosity_model(system::TotalLagrangianSPHSystem) end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, neighbor_system::TotalLagrangianSPHSystem, - density_calculator) + density_calculator, correction) (; boundary_model) = neighbor_system (; smoothing_length) = particle_system # Pressure acceleration for fluid-solid interaction. This is identical to # `pressure_acceleration` for the `BoundarySPHSystem`. - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - smoothing_length, grad_kernel, boundary_model, - density_calculator) + return pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model, density_calculator, correction) end diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index e7ee22392..18f98f117 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -6,14 +6,14 @@ integral_2d_radial, _ = quadgk(r -> r * TrixiParticles.kernel(smk, r, 1.0), 0, TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) - return 2 * π * integral_2d_radial + return 2 * pi * integral_2d_radial end function integrate_kernel_3d(smk) integral_3d_radial, _ = quadgk(r -> r^2 * TrixiParticles.kernel(smk, r, 1.0), 0, TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) - return 4 * π * integral_3d_radial + return 4 * pi * integral_3d_radial end # Treat the truncated Gaussian kernel separately diff --git a/test/schemes/fluid/weakly_compressible_sph/rhs.jl b/test/schemes/fluid/weakly_compressible_sph/rhs.jl index 6e98f1c2f..a4a6a5e8a 100644 --- a/test/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -51,13 +51,17 @@ # Compute accelerations a -> b and b -> a dv1 = TrixiParticles.pressure_acceleration(1.0, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - grad_kernel, system, system, - density_calculator) + distance, + grad_kernel, system, -1, + system, + density_calculator, nothing) dv2 = TrixiParticles.pressure_acceleration(1.0, m_a, p_b, p_a, rho_b, rho_a, -pos_diff, - -grad_kernel, system, system, - density_calculator) + distance, + -grad_kernel, system, -1, + system, + density_calculator, nothing) # Test that both forces are identical but in opposite directions @test isapprox(m_a * dv1, -m_b * dv2, rtol=2eps()) diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 99d0795bc..98ab834b6 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -88,7 +88,7 @@ "SphereShape 2D", "RectangularShape 2D with ShepardKernelCorrection", "RectangularShape 2D with AkinciFreeSurfaceCorrection", - "RectangularShape 2D with KernelGradientCorrection", + "RectangularShape 2D with KernelCorrection", ] NDIMS_ = [2, 3, 2, 3, 2, 2, 2, 2] density_calculators = [ @@ -103,7 +103,7 @@ Nothing(), ShepardKernelCorrection(), AkinciFreeSurfaceCorrection(1000.0), - KernelGradientCorrection(), + KernelCorrection(), ] @testset "$(setup_names[i])" for i in eachindex(setups) @@ -146,7 +146,7 @@ if density_calculator isa SummationDensity @test length(system.cache.density) == size(setup.coordinates, 2) end - if corr isa ShepardKernelCorrection || corr isa KernelGradientCorrection + if corr isa ShepardKernelCorrection || corr isa KernelCorrection @test length(system.cache.kernel_correction_coefficient) == size(setup.coordinates, 2) end From 15f42a6b97c63eb3f47cb1c94bc85357d7a02f6e Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 5 Jan 2024 16:54:59 +0100 Subject: [PATCH 12/12] Implement perturbed rectangular patch for better testing (#324) * Add rectangular patch * Use rectangular patch for momentum and energy tests * Fix non-perturbed center particle coordinates * Rename `initial_condition` to `fluid` * Refactor * Reformat * Add `Random` to test dependencies * Fix typo * Update test/schemes/fluid/weakly_compressible_sph/rhs.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> * Update test/schemes/fluid/weakly_compressible_sph/rhs.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> * Update test/schemes/fluid/weakly_compressible_sph/rhs.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> * Update test/schemes/fluid/weakly_compressible_sph/rhs.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --------- Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- test/Project.toml | 1 + .../fluid/weakly_compressible_sph/rhs.jl | 217 ++++++++---------- test/test_util.jl | 33 +++ 3 files changed, 131 insertions(+), 120 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 85ced8502..897dd26f0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/test/schemes/fluid/weakly_compressible_sph/rhs.jl b/test/schemes/fluid/weakly_compressible_sph/rhs.jl index a4a6a5e8a..09819e898 100644 --- a/test/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -38,9 +38,8 @@ TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 smoothing_length = -1.0 - initial_condition = InitialCondition(coordinates, velocity, mass, - density) - system = WeaklyCompressibleSPHSystem(initial_condition, + fluid = InitialCondition(coordinates, velocity, mass, density) + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, state_equation, smoothing_kernel, smoothing_length) @@ -70,142 +69,120 @@ end end - @testset verbose=true "`interact!`" begin - # The following tests for linear and angular momentum and total energy conservation - # are based on Sections 3.3.4 and 3.4.2 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 and Total Energy 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 + # The following tests for linear and angular momentum and total energy conservation + # are based on Sections 3.3.4 and 3.4.2 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 and Total Energy Conservation" begin + # We are testing the momentum conservation of SPH with random initial configurations + density_calculators = [ContinuityDensity(), SummationDensity()] + + particle_spacing = 0.1 + + # 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 = 1.2particle_spacing + search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length) + + @testset "`$(nameof(typeof(density_calculator)))`" for density_calculator in density_calculators + # Run three times with different seed for the random initial condition + for seed in 1:3 + # A larger number of particles will increase accumulated errors in the + # summation. A larger tolerance has to be used for the tests below. + fluid = rectangular_patch(particle_spacing, (3, 3), seed=seed) + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, + state_equation, smoothing_kernel, + smoothing_length) + n_particles = TrixiParticles.nparticles(system) + n_vars = TrixiParticles.v_nvariables(system) + + # Overwrite `system.pressure` because we skip the update step + system.pressure .= fluid.pressure + + u = fluid.coordinates + if density_calculator isa SummationDensity + # Density is stored in the cache + v = fluid.velocity + system.cache.density .= fluid.density + else + # Density is integrated with `ContinuityDensity` + v = vcat(fluid.velocity, fluid.density') + end - nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius, - TrixiParticles.eachparticle(system)) + nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius, + TrixiParticles.eachparticle(system)) - # Result - dv = zeros(3, 8) - TrixiParticles.interact!(dv, v, u, v, u, nhs, system, system) + # Result + dv = zeros(n_vars, n_particles) + 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) + # Linear momentum conservation + # ∑ m_a dv_a + deriv_linear_momentum = sum(fluid.mass' .* view(dv, 1:2, :), dims=2) - @test isapprox(deriv_linear_momentum, zeros(2, 1), atol=3e-14) + @test isapprox(deriv_linear_momentum, zeros(2, 1), atol=4e-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) + # 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 + return fluid.mass[particle] * cross(r_a, dv_a) + end - # ∑ m_a (r_a × dv_a) - deriv_angular_momentum = sum(deriv_angular_momentum, 1:8) + # ∑ m_a (r_a × dv_a) + deriv_angular_momentum = sum(deriv_angular_momentum, 1:n_particles) - @test isapprox(deriv_angular_momentum, zeros(3), atol=2e-14) + # Cross product is always 3-dimensional + @test isapprox(deriv_angular_momentum, zeros(3), atol=3e-15) - # Total energy conservation - drho(::ContinuityDensity, particle) = dv[3, particle] + # Total energy conservation + drho(::ContinuityDensity, particle) = dv[end, particle] - # Derivative of the density summation. This is a slightly different - # formulation of the continuity equation. - function drho_particle(particle, neighbor) - m_b = mass[i][neighbor] - vdiff = TrixiParticles.current_velocity(v, system, particle) - - TrixiParticles.current_velocity(v, system, neighbor) + function drho(::SummationDensity, particle) + return sum(neighbor -> drho_particle(particle, neighbor), 1:n_particles) + end - pos_diff = TrixiParticles.current_coords(u, system, particle) - - TrixiParticles.current_coords(u, system, neighbor) - distance = norm(pos_diff) + # Derivative of the density summation. This is a slightly different + # formulation of the continuity equation. + function drho_particle(particle, neighbor) + m_b = TrixiParticles.hydrodynamic_mass(system, neighbor) + v_diff = TrixiParticles.current_velocity(v, system, particle) - + TrixiParticles.current_velocity(v, system, neighbor) - # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return 0.0 + pos_diff = TrixiParticles.current_coords(u, system, particle) - + TrixiParticles.current_coords(u, system, neighbor) + distance = norm(pos_diff) - grad_kernel = TrixiParticles.smoothing_kernel_grad(system, - pos_diff, - distance) + # Only consider particles with a distance > 0 + distance < sqrt(eps()) && return 0.0 - return m_b * dot(vdiff, grad_kernel) - end + grad_kernel = TrixiParticles.smoothing_kernel_grad(system, pos_diff, + distance) - function drho(::SummationDensity, particle) - return sum(neighbor -> drho_particle(particle, neighbor), 1:8) - end + return m_b * dot(v_diff, grad_kernel) + end - # m_a (v_a ⋅ dv_a + dte_a), - # where `te` is the thermal energy, called `u` in the Price paper. - function deriv_energy(particle) - dte_a = pressure[i][particle] / density[i][particle]^2 * - drho(density_calculator, particle) - v_a = TrixiParticles.extract_svector(v, system, particle) - dv_a = TrixiParticles.extract_svector(dv, system, particle) + # m_a (v_a ⋅ dv_a + dte_a), + # where `te` is the thermal energy, called `u` in the Price paper. + function deriv_energy(particle) + p_a = fluid.pressure[particle] + rho_a = fluid.density[particle] + dte_a = p_a / rho_a^2 * drho(density_calculator, particle) + v_a = TrixiParticles.extract_svector(v, system, particle) + dv_a = TrixiParticles.extract_svector(dv, system, particle) - return mass[i][particle] * (dot(v_a, dv_a) + dte_a) - end + return fluid.mass[particle] * (dot(v_a, dv_a) + dte_a) + end - # ∑ m_a (v_a ⋅ dv_a + dte_a) - deriv_total_energy = sum(deriv_energy, 1:8) + # ∑ m_a (v_a ⋅ dv_a + dte_a) + deriv_total_energy = sum(deriv_energy, 1:n_particles) - @test isapprox(deriv_total_energy, 0.0, atol=4e-14) - end + @test isapprox(deriv_total_energy, 0.0, atol=2e-15) end end end diff --git a/test/test_util.jl b/test/test_util.jl index 6897bad29..c087a4dea 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -3,6 +3,7 @@ using TrixiParticles using LinearAlgebra using Printf using QuadGK: quadgk +using Random: Random """ @trixi_testset "name of the testset" #= code to test #= @@ -32,3 +33,35 @@ macro trixi_testset(name, expr) nothing end end + +function perturb!(data, amplitude) + for i in eachindex(data) + # Perturbation in the interval (-amplitude, amplitude) + data[i] += 2amplitude * rand() - amplitude + end + + return data +end + +# Rectangular patch of particles, optionally with a perturbation in position and quantities +function rectangular_patch(particle_spacing, size; density=1000.0, pressure=0.0, seed=1, + perturbation_factor=1.0) + # Fixed seed to ensure reproducibility + Random.seed!(seed) + + # Center particle at the origin (assuming odd size) + min_corner = -particle_spacing / 2 .* size + ic = RectangularShape(particle_spacing, size, min_corner, density, pressure=pressure) + + perturb!(ic.mass, perturbation_factor * 0.1 * ic.mass[1]) + perturb!(ic.density, perturbation_factor * 0.1density) + perturb!(ic.pressure, perturbation_factor * 2000) + perturb!(ic.velocity, perturbation_factor * 0.5particle_spacing) + perturb!(ic.coordinates, perturbation_factor * 0.5particle_spacing) + + # Don't perturb center particle position + center_particle = ceil(Int, prod(size) / 2) + ic.coordinates[:, center_particle] .= 0.0 + + return ic +end