Merge branch 'main' into validate-open-boundaries

Author: svchb

Project.toml

@@ -50,7 +50,7 @@ FastPow = "0.1"
 FileIO = "1"
 ForwardDiff = "1"
 GPUArraysCore = "0.2"
-JSON = "0.21"
+JSON = "1"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
 OrdinaryDiffEq = "6.91"

src/TrixiParticles.jl

@@ -15,7 +15,7 @@ using ForwardDiff: ForwardDiff
 using GPUArraysCore: AbstractGPUArray
 using JSON: JSON
 using KernelAbstractions: KernelAbstractions, @kernel, @index
-using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
+using LinearAlgebra: norm, normalize, cross, dot, I, tr, inv, pinv, det
 using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf

src/general/interpolation.jl

@@ -596,6 +596,12 @@ end
                                  distance, grad_kernel, correction)
 end
 
+# The function above misuses the pressure acceleration function by passing a Matrix as `p_a`.
+# This doesn't work with `tensile_instability_control`, so we disable TIC in this case.
+@inline function tensile_instability_control(m_a, m_b, rho_a, rho_b, p_a::SMatrix, p_b, W_a)
+    return pressure_acceleration_continuity_density(m_a, m_b, rho_a, rho_b, p_a, A_b, W_a)
+end
+
 function continuity_equation_shifting!(dv, shifting::TransportVelocityAdami{true},
                                        particle_system, neighbor_system,
                                        particle, neighbor, grad_kernel, rho_a, rho_b, m_b)

src/schemes/fluid/shifting_techniques.jl

@@ -17,11 +17,11 @@ TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26"
 
 [compat]
 CSV = "0.10"
-CairoMakie = "0.13, 0.15"
+CairoMakie = "0.15"
 DataFrames = "1.6"
 GLM = "1.9"
 Glob = "1.3"
-JSON = "0.21"
+JSON = "1"
 OrdinaryDiffEq = "6.49"
 Plots = "1"
 Polyester = "0.7"

test/Project.toml

@@ -107,6 +107,7 @@
 
         boundary_model = BoundaryModelDummyParticles(bnd.density, bnd.mass,
                                                      state_equation=state_equation,
+                                                     viscosity=viscosity,
                                                      AdamiPressureExtrapolation(),
                                                      smoothing_kernel, smoothing_length)
 
@@ -552,6 +553,42 @@
                 end
             end
         end
+
+        @testset verbose=true "Include Wall Velocity" begin
+            # Define vertical interpolation line
+            start_svector = SVector(0.0, 0.0)
+            end_svector = SVector(0.0, 1.0)
+            points_coords_ = range(start_svector, end_svector, length=10)
+            points_coords = collect(reinterpret(reshape, eltype(start_svector),
+                                                points_coords_))
+
+            # Linear velocity field for fluid and boundary
+            v_fluid = InitialCondition(;
+                                       coordinates=fluid_system.initial_condition.coordinates,
+                                       density=1000.0, particle_spacing,
+                                       velocity=(pos) -> SVector(0.0, pos[2])).velocity
+            v_boundary = InitialCondition(; coordinates=boundary_system.coordinates,
+                                          density=1000.0, particle_spacing,
+                                          velocity=(pos) -> SVector(0.0, pos[2])).velocity
+
+            boundary_system.boundary_model.cache.wall_velocity .= v_boundary
+
+            v_ode = vcat(v_fluid, fluid.density')
+            u_ode = fluid_system.initial_condition.coordinates
+
+            v_wall_velocity = interpolate_points(points_coords, semi_boundary,
+                                                 include_wall_velocity=true,
+                                                 fluid_system, v_ode, u_ode).velocity
+
+            v_no_wall_velocity = interpolate_points(points_coords, semi_boundary,
+                                                    include_wall_velocity=false,
+                                                    fluid_system, v_ode, u_ode).velocity
+
+            @test isapprox(v_wall_velocity[2, 1], 0.0; atol=eps())
+            @test isapprox(v_no_wall_velocity[2, 1], 0.1; atol=eps())
+            @test any(isapprox.(v_wall_velocity[:, 3:end], v_no_wall_velocity[:, 3:end],
+                                atol=eps()))
+        end
     end
 
     @testset verbose=true "3D" begin
@@ -621,6 +658,7 @@
 
         boundary_model = BoundaryModelDummyParticles(bnd.density, bnd.mass,
                                                      state_equation=state_equation,
+                                                     viscosity=viscosity,
                                                      AdamiPressureExtrapolation(),
                                                      smoothing_kernel, smoothing_length)
 
@@ -883,5 +921,41 @@
                 compare_interpolation_result(result, expected_res)
             end
         end
+
+        @testset verbose=true "Include Wall Velocity" begin
+            # Define vertical interpolation line
+            start_svector = SVector(0.0, 0.0, 0.0)
+            end_svector = SVector(0.0, 1.0, 0.0)
+            points_coords_ = range(start_svector, end_svector, length=10)
+            points_coords = collect(reinterpret(reshape, eltype(start_svector),
+                                                points_coords_))
+
+            # Linear velocity field for fluid and boundary
+            v_fluid = InitialCondition(;
+                                       coordinates=fluid_system.initial_condition.coordinates,
+                                       density=1000.0, particle_spacing,
+                                       velocity=(pos) -> SVector(0.0, pos[2], 0.0)).velocity
+            v_boundary = InitialCondition(; coordinates=boundary_system.coordinates,
+                                          density=1000.0, particle_spacing,
+                                          velocity=(pos) -> SVector(0.0, pos[2], 0.0)).velocity
+
+            boundary_system.boundary_model.cache.wall_velocity .= v_boundary
+
+            v_ode = vcat(v_fluid, fluid.density')
+            u_ode = fluid_system.initial_condition.coordinates
+
+            v_wall_velocity = interpolate_points(points_coords, semi_boundary,
+                                                 include_wall_velocity=true,
+                                                 fluid_system, v_ode, u_ode).velocity
+
+            v_no_wall_velocity = interpolate_points(points_coords, semi_boundary,
+                                                    include_wall_velocity=false,
+                                                    fluid_system, v_ode, u_ode).velocity
+
+            @test isapprox(v_wall_velocity[2, 1], 0.0; atol=eps())
+            @test isapprox(v_no_wall_velocity[2, 1], 0.1; atol=eps())
+            @test any(isapprox.(v_wall_velocity[:, 3:end], v_no_wall_velocity[:, 3:end],
+                                atol=eps()))
+        end
     end
 end;

test/general/interpolation.jl

@@ -253,7 +253,7 @@
                                                  time_step=0.001)
 
         u0 = zeros(TrixiParticles.u_nvariables(system),
-                   TrixiParticles.n_moving_particles(system))
+                   TrixiParticles.n_integrated_particles(system))
         TrixiParticles.write_u0!(u0, system)
 
         @test u0 == coordinates
@@ -279,7 +279,7 @@
                                                  time_step=0.001)
 
         v0 = zeros(TrixiParticles.v_nvariables(system),
-                   TrixiParticles.n_moving_particles(system))
+                   TrixiParticles.n_integrated_particles(system))
         TrixiParticles.write_v0!(v0, system)
 
         @test v0 == velocity

test/systems/iisph_system.jl

@@ -5,3 +5,4 @@ include("boundary_system.jl")
 include("open_boundary_system.jl")
 include("dem_system.jl")
 include("packing_system.jl")
+include("iisph_system.jl")