Merge remote-tracking branch 'origin/generalize_surface_normal_calc' …

…into morris_surface_tension

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/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="")
 
 # ==========================================================================================

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/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

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)

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/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")