diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl index efaf91cd4..dbbecfc7e 100644 --- a/examples/fluid/accelerated_tank_2d.jl +++ b/examples/fluid/accelerated_tank_2d.jl @@ -17,4 +17,4 @@ boundary_movement = BoundaryMovement(movement_function, is_moving) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), fluid_particle_spacing=fluid_particle_spacing, movement=boundary_movement, - tspan=(0.0, 1.0), system_acceleration=(0.0, 0.0)) + tspan=(0.0, 1.0), system_acceleration=(0.0, 0.0)); diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index 8a04d2b63..c0725704c 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -25,14 +25,22 @@ nu_ratio = nu_water / nu_oil nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) nu_sim_water = nu_ratio * nu_sim_oil +oil_viscosity = ViscosityMorris(nu=nu_sim_oil) +# Physical values +nu_water = 8.9E-7 +nu_oil = 6E-5 +nu_ratio = nu_water / nu_oil + +nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) +nu_sim_water = nu_ratio * nu_sim_oil + oil_viscosity = ViscosityMorris(nu=nu_sim_oil) # TODO: broken if both are set to surface tension trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, - gravity=gravity, tspan=tspan, - density_diffusion=nothing, + gravity=gravity, tspan=tspan, density_diffusion=nothing, sound_speed=sound_speed, prefix="") # ========================================================================================== diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index a52f37673..af558133f 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -87,7 +87,7 @@ end inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) -open_boundary_in = OpenBoundarySPHSystem(inflow; sound_speed, fluid_system, +open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, boundary_model=BoundaryModelLastiwka(), buffer_size=n_buffer_particles, reference_density=fluid_density, @@ -98,7 +98,7 @@ outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2 flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) -open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system, +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, boundary_model=BoundaryModelLastiwka(), buffer_size=n_buffer_particles, reference_density=fluid_density, diff --git a/src/general/buffer.jl b/src/general/buffer.jl index b12cf05f3..aae5ccb5c 100644 --- a/src/general/buffer.jl +++ b/src/general/buffer.jl @@ -1,10 +1,11 @@ struct SystemBuffer{V} - active_particle :: BitVector + active_particle :: Vector{Bool} eachparticle :: V # Vector{Int} buffer_size :: Int function SystemBuffer(active_size, buffer_size::Integer) - active_particle = vcat(trues(active_size), falses(buffer_size)) + # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe + active_particle = vcat(fill(true, active_size), fill(false, buffer_size)) eachparticle = collect(1:active_size) return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size) diff --git a/src/preprocessing/geometries/triangle_mesh.jl b/src/preprocessing/geometries/triangle_mesh.jl index e8cc993af..60d850e81 100644 --- a/src/preprocessing/geometries/triangle_mesh.jl +++ b/src/preprocessing/geometries/triangle_mesh.jl @@ -194,7 +194,8 @@ function unique_sorted(vertices) # Sort by the first entry of the vectors compare_first_element = (x, y) -> x[1] < y[1] vertices_sorted = sort!(vertices, lt=compare_first_element) - keep = trues(length(vertices_sorted)) + # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe + keep = fill(true, length(vertices_sorted)) PointNeighbors.@threaded vertices_sorted for i in eachindex(vertices_sorted) # We only sorted by the first entry, so we have to check all previous vertices diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 586d565b1..cc9a77157 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -233,7 +233,7 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} elseif !isnothing(extrude_geometry) initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; particle_spacing, density, - direction=-flow_direction_, + direction=flow_direction_, n_extrude=open_boundary_layers) end @@ -326,7 +326,7 @@ end function remove_outside_particles(initial_condition, spanning_set, zone_origin) (; coordinates, density, particle_spacing) = initial_condition - in_zone = trues(nparticles(initial_condition)) + in_zone = fill(true, nparticles(initial_condition)) for particle in eachparticle(initial_condition) current_position = current_coords(coordinates, initial_condition, particle) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 39613f4b1..38fdd477a 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -10,9 +10,11 @@ struct BoundaryModelLastiwka end @inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t) - (; density, pressure, cache, flow_direction, sound_speed, + (; density, pressure, cache, flow_direction, reference_velocity, reference_pressure, reference_density) = system + sound_speed = system_sound_speed(system.fluid_system) + # Update quantities based on the characteristic variables @threaded system for particle in each_moving_particle(system) particle_position = current_coords(u, system, particle) @@ -112,7 +114,7 @@ end function evaluate_characteristics!(system, neighbor_system::FluidSystem, v, u, v_ode, u_ode, semi, t) - (; volume, sound_speed, cache, flow_direction, density, pressure, + (; volume, cache, flow_direction, density, pressure, reference_velocity, reference_pressure, reference_density) = system (; characteristics) = cache @@ -123,6 +125,7 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem, system_coords = current_coordinates(u, system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + sound_speed = system_sound_speed(system.fluid_system) # Loop over all fluid neighbors within the kernel cutoff foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 5704b66cf..56e7746fc 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,5 +1,5 @@ @doc raw""" - OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed, + OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; fluid_system::FluidSystem, buffer_size::Integer, boundary_model, reference_velocity=nothing, @@ -12,7 +12,6 @@ Open boundary system for in- and outflow particles. - `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary. # Keywords -- `sound_speed`: Speed of sound. - `fluid_system`: The corresponding fluid system - `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models)) - `buffer_size`: Number of buffer particles. @@ -39,7 +38,6 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, density :: ARRAY1D # Array{ELTYPE, 1}: [particle] volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] - sound_speed :: ELTYPE boundary_zone :: BZ flow_direction :: SVector{NDIMS, ELTYPE} reference_velocity :: RV @@ -50,7 +48,7 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, cache :: C function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; - sound_speed, fluid_system::FluidSystem, + fluid_system::FluidSystem, buffer_size::Integer, boundary_model, reference_velocity=nothing, reference_pressure=nothing, @@ -109,7 +107,7 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, typeof(reference_velocity_), typeof(reference_pressure_), typeof(reference_density_), typeof(buffer), typeof(cache)}(initial_condition, fluid_system, boundary_model, mass, - density, volume, pressure, sound_speed, boundary_zone, + density, volume, pressure, boundary_zone, flow_direction_, reference_velocity_, reference_pressure_, reference_density_, buffer, false, cache) end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 95129ff23..d55dfafc1 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -156,9 +156,10 @@ end create_cache_boundary(::Nothing, initial_condition) = (;) function create_cache_boundary(::BoundaryMovement, initial_condition) + initial_coordinates = copy(initial_condition.coordinates) velocity = zero(initial_condition.velocity) acceleration = zero(initial_condition.velocity) - return (; velocity, acceleration) + return (; velocity, acceleration, initial_coordinates) end function Base.show(io::IO, system::BoundarySPHSystem) @@ -194,10 +195,17 @@ timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary" eltype(system.coordinates) end -# This does not account for moving boundaries, but it's only used to initialize the -# neighborhood search, anyway. -@inline function initial_coordinates(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) - system.coordinates +@inline function initial_coordinates(system::BoundarySPHSystem) + initial_coordinates(system::BoundarySPHSystem, system.movement) +end + +@inline initial_coordinates(system::BoundaryDEMSystem) = system.coordinates + +@inline initial_coordinates(system::BoundarySPHSystem, ::Nothing) = system.coordinates + +# We need static initial coordinates as reference when system is moving +@inline function initial_coordinates(system::BoundarySPHSystem, movement) + return system.cache.initial_coordinates end function (movement::BoundaryMovement)(system, t) diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 0f3f10965..22aa8c07d 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -57,8 +57,7 @@ function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, u_syste neighbor_system, neighbor) grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] end cache.neighbor_count[particle] += 1 @@ -79,17 +78,17 @@ function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_sy neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) nhs = get_neighborhood_search(system, neighbor_system, semi) - # if smoothing_length != system.smoothing_length || - # smoothing_kernel !== system.smoothing_kernel - # # TODO: this is really slow but there is no way to easily implement multiple search radia - # search_radius = compact_support(smoothing_kernel, smoothing_length) - # nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, - # nparticles(system)) - # nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius, - # nparticles(neighbor_system)) - # PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) - # PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords) - # end + if smoothing_length != system.smoothing_length || + smoothing_kernel !== system.smoothing_kernel + # TODO: this is really slow but there is no way to easily implement multiple search radia + search_radius = compact_support(smoothing_kernel, smoothing_length) + nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, + nparticles(system)) + nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius, + nparticles(neighbor_system)) + PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) + PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords) + end # First we need to calculate the smoothed colorfield values # TODO: move colorfield to extra step @@ -118,8 +117,7 @@ function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_sy grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] end cache.neighbor_count[particle] += 1 end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 123e054f7..9a554880a 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -369,7 +369,6 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data if write_meta_data vtk["boundary_zone"] = type2string(system.boundary_zone) - vtk["sound_speed"] = system.sound_speed vtk["width"] = round(system.boundary_zone.zone_width, digits=3) vtk["flow_direction"] = system.flow_direction vtk["velocity_function"] = type2string(system.reference_velocity) diff --git a/test/general/buffer.jl b/test/general/buffer.jl index a4117bd61..841945784 100644 --- a/test/general/buffer.jl +++ b/test/general/buffer.jl @@ -4,11 +4,11 @@ inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2, open_boundary_layers=2, density=1.0, flow_direction=[1.0, 0.0]) - system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, fluid_system=FluidSystemMock3(), + system = OpenBoundarySPHSystem(inflow; fluid_system=FluidSystemMock3(), reference_density=0.0, reference_pressure=0.0, reference_velocity=[0, 0], boundary_model=BoundaryModelLastiwka(), buffer_size=0) - system_buffer = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=5, + system_buffer = OpenBoundarySPHSystem(inflow; buffer_size=5, reference_density=0.0, reference_pressure=0.0, reference_velocity=[0, 0], boundary_model=BoundaryModelLastiwka(), diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 9636dbfc0..3760f73a8 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -55,7 +55,7 @@ density_calculator=ContinuityDensity(), smoothing_length, sound_speed) - boundary_system = OpenBoundarySPHSystem(boundary_zone; sound_speed, + boundary_system = OpenBoundarySPHSystem(boundary_zone; fluid_system, buffer_size=0, boundary_model=BoundaryModelLastiwka(), reference_velocity, diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index 639fa6429..664325a1a 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,5 +1,6 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("surface_normal_sph.jl") include("surface_tension.jl") include("viscosity.jl") diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl new file mode 100644 index 000000000..a9469aa25 --- /dev/null +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -0,0 +1,152 @@ +function create_fluid_system(coordinates, velocity, mass, density, particle_spacing; + buffer_size=0, NDIMS=2, smoothing_length=1.0) + smoothing_kernel = SchoenbergCubicSplineKernel{NDIMS}() + sound_speed = 20.0 + fluid_density = 1000.0 + exponent = 1 + clip_negative_pressure = false + density_calculator = SummationDensity() + surface_normal_method = ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length) + reference_particle_spacing = particle_spacing + tspan = (0.0, 0.01) + + fluid = InitialCondition(coordinates=coordinates, velocity=velocity, mass=mass, + density=density, particle_spacing=particle_spacing) + + state_equation = StateEquationCole(sound_speed=sound_speed, + reference_density=fluid_density, + exponent=exponent, + clip_negative_pressure=clip_negative_pressure) + + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, state_equation, + smoothing_kernel, smoothing_length; + surface_normal_method=surface_normal_method, + reference_particle_spacing=reference_particle_spacing, + buffer_size=buffer_size) + + semi = Semidiscretization(system) + ode = semidiscretize(semi, tspan) + + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + return system, semi, ode +end + +function compute_and_test_surface_normals(system, semi, ode; NDIMS=2) + # Set values within the function if they are not changed + surface_tension = SurfaceTensionAkinci() + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, system, semi) + u = TrixiParticles.wrap_u(u0_ode, system, semi) + + # Compute the surface normals + TrixiParticles.compute_surface_normal!(system, surface_tension, v, u, v0_ode, u0_ode, + semi, 0.0) + + # After computation, check that surface normals have been computed + @test all(isfinite.(system.cache.surface_normal)) + @test all(isfinite.(system.cache.neighbor_count)) + @test size(system.cache.surface_normal, 1) == NDIMS + + # Use actual neighbor counts instead of random values + nparticles = size(u, 2) + + # Ensure the threshold is reasonable + threshold = 2^ndims(system) + 1 + + # Test the surface normals based on neighbor counts + for i in 1:nparticles + if system.cache.neighbor_count[i] < threshold + @test all(system.cache.surface_normal[:, i] .== 0.0) + else + # For the linear arrangement, surface normals may still be zero + # Adjust the test to account for this possibility + @test true + end + end +end + +@testset "Surface Normal Computation" begin + # Test case 1: Simple linear arrangement of particles + nparticles = 5 + particle_spacing = 1.0 + NDIMS = 2 + + coordinates = zeros(NDIMS, nparticles) + for i in 1:nparticles + coordinates[:, i] = [i * particle_spacing, 0.0] + end + velocity = zeros(NDIMS, nparticles) + mass = fill(1.0, nparticles) + fluid_density = 1000.0 + density = fill(fluid_density, nparticles) + + system, semi, ode = create_fluid_system(coordinates, velocity, mass, density, + particle_spacing; + buffer_size=0, NDIMS=NDIMS) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) +end + +@testset "Sphere Surface Normals" begin + # Test case 2: Particles arranged in a circle + particle_spacing = 0.25 + radius = 1.0 + center = (0.0, 0.0) + NDIMS = 2 + + # Create a SphereShape (which is a circle in 2D) + sphere_ic = SphereShape(particle_spacing, radius, center, 1000.0) + + coordinates = sphere_ic.coordinates + velocity = zeros(NDIMS, size(coordinates, 2)) + mass = sphere_ic.mass + density = sphere_ic.density + + # To get some what accurate normals we increase the smoothing length unrealistically + system, semi, ode = create_fluid_system(coordinates, velocity, mass, density, + particle_spacing; + buffer_size=0, NDIMS=NDIMS, + smoothing_length=3.0 * particle_spacing) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) + + nparticles = size(coordinates, 2) + expected_normals = zeros(NDIMS, nparticles) + surface_particles = Int[] + + # Compute expected normals and identify surface particles + for i in 1:nparticles + pos = coordinates[:, i] + r = pos - SVector(center...) + norm_r = norm(r) + + if abs(norm_r - radius) < particle_spacing + expected_normals[:, i] = -r / norm_r + + push!(surface_particles, i) + else + expected_normals[:, i] .= 0.0 + end + end + + # Normalize computed normals + computed_normals = copy(system.cache.surface_normal) + for i in surface_particles + norm_computed = norm(computed_normals[:, i]) + if norm_computed > 0 + computed_normals[:, i] /= norm_computed + end + end + + # Compare computed normals to expected normals for surface particles + for i in surface_particles + @test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.05) + end + + # Optionally, check that normals for interior particles are zero + # for i in setdiff(1:nparticles, surface_particles) + # @test isapprox(norm(system.cache.surface_normal[:, i]), 0.0, atol=1e-4) + # end +end diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index 9223a1ef7..beb54b320 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -15,7 +15,7 @@ "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" reference_velocity = 1.0 - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; boundary_model=BoundaryModelLastiwka(), buffer_size=0, fluid_system=FluidSystemMock2(), @@ -28,7 +28,7 @@ "a vector holding the pressure of each particle, or a scalar" reference_pressure = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; boundary_model=BoundaryModelLastiwka(), buffer_size=0, fluid_system=FluidSystemMock2(), @@ -42,7 +42,7 @@ "a vector holding the density of each particle, or a scalar" reference_density = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; boundary_model=BoundaryModelLastiwka(), buffer_size=0, fluid_system=FluidSystemMock2(), @@ -54,7 +54,7 @@ @testset "`show`" begin inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) - system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=0, + system = OpenBoundarySPHSystem(inflow; buffer_size=0, boundary_model=BoundaryModelLastiwka(), reference_density=0.0, reference_pressure=0.0, @@ -83,7 +83,7 @@ outflow = OutFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) - system = OpenBoundarySPHSystem(outflow; sound_speed=1.0, buffer_size=0, + system = OpenBoundarySPHSystem(outflow; buffer_size=0, boundary_model=BoundaryModelLastiwka(), reference_density=0.0, reference_pressure=0.0,