diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index af558133f..c8778f485 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -48,7 +48,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi # Shift pipe walls in negative x-direction for the inflow pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers -n_buffer_particles = 4 * pipe.n_particles_per_dimension[2] +n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1) # ========================================================================================== # ==== Fluid @@ -77,14 +77,16 @@ 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) return SVector(prescribed_velocity, 0.0) end -inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, +plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = InFlow(; plane=plane_in, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, @@ -92,9 +94,10 @@ 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]]), +plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) +outflow = OutFlow(; plane=plane_out, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) @@ -103,7 +106,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 diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl new file mode 100644 index 000000000..fb2cc7f0a --- /dev/null +++ b/examples/fluid/pipe_flow_3d.jl @@ -0,0 +1,45 @@ +# 3D channel flow simulation with open boundaries. + +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +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]) + +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, boundary_size=boundary_size, + flow_direction=flow_direction, faces=(false, false, true, true, true, true), + tspan=tspan, smoothing_kernel=WendlandC2Kernel{3}(), + reference_velocity=velocity_function3d, + plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0], + [0.0, 0.0, domain_size[3]]), + plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0], + [domain_size[1], 0.0, domain_size[3]])) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index cc9a77157..ad72045d6 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/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 diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 38fdd477a..f803bf9ca 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/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 + characteristics[2, particle] = 0 + characteristics[3, particle] = 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 diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9d2c79b00..c8d80ca8f 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/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) diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index c666ec3d7..9b31c463c 100644 --- a/src/setups/extrude_geometry.jl +++ b/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 diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 2e5462a35..3673e5bd5 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -238,6 +238,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", diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index ad1b69c4b..2ed6e2810 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/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,