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/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/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index fbc5bb336..9b33406a2 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -336,7 +336,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/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,