Fixes for Open Boundary Systems for free surfaces

examples/fluid/pipe_flow_2d.jl

@@ -77,7 +77,7 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi
 
 # ==========================================================================================
 # ==== Open Boundary
-function velocity_function(pos, t)
+function velocity_function2d(pos, t)
  # Use this for a time-dependent inflow velocity
  # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0)
 
@@ -92,7 +92,7 @@ open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
  buffer_size=n_buffer_particles,
  reference_density=fluid_density,
  reference_pressure=pressure,
- reference_velocity=velocity_function)
+ reference_velocity=velocity_function2d)
 
 outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
  flow_direction, open_boundary_layers, density=fluid_density,
@@ -103,7 +103,7 @@ open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
  buffer_size=n_buffer_particles,
  reference_density=fluid_density,
  reference_pressure=pressure,
- reference_velocity=velocity_function)
+ reference_velocity=velocity_function2d)
 
 # ==========================================================================================
 # ==== Boundary

examples/fluid/pipe_flow_3d.jl

@@ -0,0 +1,46 @@
+# 3D channel flow simulation with open boundaries.
+using TrixiParticles
+using OrdinaryDiffEq
+
+tspan = (0.0, 2.0)
+
+# load variables
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"),
+ sol=nothing)
+
+function velocity_function3d(pos, t)
+ # Use this for a time-dependent inflow velocity
+ # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0)
+
+ return SVector(prescribed_velocity, 0.0, 0.0)
+end
+
+domain_size = (1.0, 0.4, 0.4)
+
+boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+ domain_size[2], domain_size[3])
+
+pipe3d = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
+ pressure=pressure, n_layers=boundary_layers,
+ faces=(false, false, true, true, true, true))
+
+flow_direction = [1.0, 0.0, 0.0]
+
+# setup simulation
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"),
+ domain_size=domain_size, flow_direction=flow_direction,
+ pipe=pipe3d,
+ n_buffer_particles=4 * pipe3d.n_particles_per_dimension[2]^2,
+ smoothing_kernel=WendlandC2Kernel{3}(),
+ reference_velocity=velocity_function3d,
+ inflow=InFlow(;
+ plane=([0.0, 0.0, 0.0],
+ [0.0, domain_size[2], 0.0],
+ [0.0, 0.0, domain_size[3]]), flow_direction,
+ open_boundary_layers, density=fluid_density, particle_spacing),
+ outflow=OutFlow(;
+ plane=([domain_size[1], 0.0, 0.0],
+ [domain_size[1], domain_size[2], 0.0],
+ [domain_size[1], 0.0, domain_size[3]]),
+ flow_direction, open_boundary_layers, density=fluid_density,
+ particle_spacing))

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -116,11 +116,11 @@ struct InFlow{NDIMS, IC, S, ZO, ZW, FD}
 
  if !isapprox(abs(dot_), 1.0, atol=1e-7)
  throw(ArgumentError("`flow_direction` is not normal to inflow plane"))
- else
- # Flip the normal vector to point in the opposite direction of `flow_direction`
- spanning_set[:, 1] .*= -sign(dot_)
  end
 
+ # Flip the normal vector to point in the opposite direction of `flow_direction`
+ spanning_set[:, 1] .*= -sign(dot_)
+
  spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set)
 
  # Remove particles outside the boundary zone.
@@ -251,11 +251,11 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD}
 
  if !isapprox(abs(dot_), 1.0, atol=1e-7)
  throw(ArgumentError("`flow_direction` is not normal to outflow plane"))
- else
- # Flip the normal vector to point in `flow_direction`
- spanning_set[:, 1] .*= sign(dot_)
  end
 
+ # Flip the normal vector to point in `flow_direction`
+ spanning_set[:, 1] .*= sign(dot_)
+
  spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set)
 
  # Remove particles outside the boundary zone.
@@ -289,8 +289,9 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width)
  edge1 = plane_points[2] - plane_points[1]
  edge2 = plane_points[3] - plane_points[1]
 
- if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7)
- throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal"))
+ # Check if the edges are linearly dependent (to avoid degenerate planes)
+ if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps())
+ throw(ArgumentError("The vectors `AB` and `AC` must not be collinear"))
  end
 
  # Calculate normal vector of plane

src/schemes/boundary/open_boundary/method_of_characteristics.jl

@@ -8,8 +8,9 @@ about the method see [description below](@ref method_of_characteristics).
 """
 struct BoundaryModelLastiwka end
 
-@inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka,
- v, u, v_ode, u_ode, semi, t)
+# Called from update callback via update_open_boundary_eachstep!
+@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka,
+ v, u, v_ode, u_ode, semi, t)
  (; density, pressure, cache, flow_direction,
  reference_velocity, reference_pressure, reference_density) = system
 
@@ -45,10 +46,13 @@ struct BoundaryModelLastiwka end
  return system
 end
 
+# Called from semidiscretization
 function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t)
  @trixi_timeit timer() "evaluate characteristics" begin
  evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
  end
+
+ return system
 end
 
 # ==== Characteristics
@@ -98,9 +102,16 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
  end
  end
 
- characteristics[1, particle] = avg_J1 / counter
- characteristics[2, particle] = avg_J2 / counter
- characteristics[3, particle] = avg_J3 / counter
+ # To prevent NANs here if the boundary has not been in contact before.
+ if counter > 0
+ characteristics[1, particle] = avg_J1 / counter
+ characteristics[2, particle] = avg_J2 / counter
+ characteristics[3, particle] = avg_J3 / counter
+ else
+ characteristics[1, particle] = 0.0
+ characteristics[2, particle] = 0.0
+ characteristics[3, particle] = 0.0
+ end
  else
  characteristics[1, particle] /= volume[particle]
  characteristics[2, particle] /= volume[particle]
@@ -164,7 +175,7 @@ end
 
 @inline function prescribe_conditions!(characteristics, particle, ::OutFlow)
  # J3 is prescribed (i.e. determined from the exterior of the domain).
- # J1 and J2 is transimtted from the domain interior.
+ # J1 and J2 is transmitted from the domain interior.
  characteristics[3, particle] = zero(eltype(characteristics))
 
  return characteristics

src/schemes/boundary/open_boundary/system.jl

@@ -77,6 +77,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  "an array where the ``i``-th column holds the velocity of particle ``i`` " *
  "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity"))
  else
+ if reference_velocity isa Function
+ test_result = reference_velocity(zeros(NDIMS), 0.0)
+ if length(test_result) != NDIMS
+ throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS"))
+ end
+ end
  reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS))
  end
 
@@ -86,6 +92,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  "each particle's coordinates and time to its pressure, " *
  "a vector holding the pressure of each particle, or a scalar"))
  else
+ if reference_pressure isa Function
+ test_result = reference_pressure(zeros(NDIMS), 0.0)
+ if length(test_result) != 1
+ throw(ArgumentError("`reference_pressure` function must be a scalar function"))
+ end
+ end
  reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS))
  end
 
@@ -95,6 +107,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV,
  "each particle's coordinates and time to its density, " *
  "a vector holding the density of each particle, or a scalar"))
  else
+ if reference_density isa Function
+ test_result = reference_density(zeros(NDIMS), 0.0)
+ if length(test_result) != 1
+ throw(ArgumentError("`reference_density` function must be a scalar function"))
+ end
+ end
  reference_density_ = wrap_reference_function(reference_density, Val(NDIMS))
  end
 
@@ -190,10 +208,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_
 
  # Update density, pressure and velocity based on the characteristic variables.
  # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971
- @trixi_timeit timer() "update quantities" update_quantities!(system,
- system.boundary_model,
- v, u, v_ode, u_ode,
- semi, t)
+ @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system,
+ system.boundary_model,
+ v, u,
+ v_ode,
+ u_ode,
+ semi, t)
 
  @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
 

src/setups/extrude_geometry.jl

@@ -177,8 +177,8 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing)
  edge2 = point3_ - point1_
 
  # Check if the points are collinear
- if norm(cross(edge1, edge2)) == 0
- throw(ArgumentError("the points must not be collinear"))
+ if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps())
+ throw(ArgumentError("The vectors `AB` and `AC` must not be collinear"))
  end
 
  # Determine the number of points along each edge

test/examples/examples.jl

@@ -226,6 +226,14 @@
  @test count_rhs_allocations(sol, semi) == 0
  end
 
+ @trixi_testset "fluid/pipe_flow_3d.jl" begin
+ @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5),
+ joinpath(examples_dir(), "fluid",
+ "pipe_flow_3d.jl"))
+ @test sol.retcode == ReturnCode.Success
+ @test count_rhs_allocations(sol, semi) == 0
+ end
+
  @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin
  @test_nowarn_mod trixi_include(@__MODULE__,
  joinpath(examples_dir(), "fluid",

test/schemes/boundary/open_boundary/boundary_zone.jl

@@ -187,10 +187,10 @@
  end
 
  @testset verbose=true "Illegal Inputs" begin
- no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]]
+ no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.4, 0.9, -0.15]]
  flow_direction = [0.0, 0.0, 1.0]
 
- error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal"
+ error_str = "The vectors `AB` and `AC` must not be collinear"
 
  @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane,
  particle_spacing=0.1,