Merge branch 'main' into generalize_surface_normal_calc

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));

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,

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)

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

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)

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,

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

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)

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)

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(),

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,

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,