Fixes for Open Boundary Systems for free surfaces

examples/fluid/pipe_flow_2d.jl

@@ -77,11 +77,15 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi
 
 # ==========================================================================================
 # ==== Open Boundary
-function velocity_function(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)
+if !@isdefined(VELOCITY_FUNC)
+    # to avoid double declaration with `trixi_include`
+    function velocity_function2d(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)
+    end
+    const VELOCITY_FUNC = true
 end
 
 inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
@@ -92,7 +96,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 +107,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/lid_driven_cavity_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,