From ed83ea1d6312476c8bc57060d5e149c36e7ccbef Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 8 Aug 2024 11:36:45 +0200 Subject: [PATCH 01/34] set test up for 1.11 --- .github/workflows/ci.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b115c6937..04070512f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -43,6 +43,7 @@ jobs: matrix: version: - '1.9' + - '1.10' - '1' os: - ubuntu-latest @@ -50,10 +51,11 @@ jobs: - x64 include: # Also run tests on Windows and macOS-ARM, but only with the latest Julia version - - version: '1' + # TODO: change back to '1' after JuliaLang/julia#55009 has been backported + - version: '1.10' os: windows-latest arch: x64 - - version: '1' + - version: '1.10' os: macos-14 arch: arm64 From bf4bea248a59bd3902131aa8061b0b71a334aab9 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 Aug 2024 12:30:59 +0200 Subject: [PATCH 02/34] fix --- examples/fluid/accelerated_tank_2d.jl | 2 +- src/schemes/boundary/system.jl | 18 +++++++++++++----- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl index efaf91cd4..dbbecfc7e 100644 --- a/examples/fluid/accelerated_tank_2d.jl +++ b/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)); diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 8f181ce7e..c097e45ab 100644 --- a/src/schemes/boundary/system.jl +++ b/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,12 +195,19 @@ 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::BoundarySPHSystem, nothing) = sytem.coordinates + +@inline function initial_coordinates(system::BoundarySPHSystem, movement) + + return system.cache.initial_coordinates +end + +@inline initial_coordinates(system::BoundaryDEMSystem) = system.coordinates + function (movement::BoundaryMovement)(system, t) (; coordinates, cache) = system (; movement_function, is_moving, moving_particles) = movement From f62923972bb0b6a2b6d58545becb5b275278cba0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 Aug 2024 12:32:48 +0200 Subject: [PATCH 03/34] typo --- src/schemes/boundary/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index c097e45ab..f057e0ebe 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -199,7 +199,7 @@ end initial_coordinates(system::BoundarySPHSystem, system.movement) end -@inline initial_coordinates(system::BoundarySPHSystem, nothing) = sytem.coordinates +@inline initial_coordinates(system::BoundarySPHSystem, nothing) = system.coordinates @inline function initial_coordinates(system::BoundarySPHSystem, movement) From 7d5a5fc2dbd95ebc5e28d532f31806b224d889b5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 Aug 2024 12:41:04 +0200 Subject: [PATCH 04/34] fix again --- src/schemes/boundary/system.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index f057e0ebe..583f3ec3e 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -199,10 +199,9 @@ end initial_coordinates(system::BoundarySPHSystem, system.movement) end -@inline initial_coordinates(system::BoundarySPHSystem, nothing) = system.coordinates +@inline initial_coordinates(system::BoundarySPHSystem, ::Nothing) = system.coordinates @inline function initial_coordinates(system::BoundarySPHSystem, movement) - return system.cache.initial_coordinates end From 9f7d1c1579d1269abd21fac0546bb0b2c92ef063 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Aug 2024 17:29:41 +0200 Subject: [PATCH 05/34] remove soundspeed from OBS --- examples/fluid/pipe_flow_2d.jl | 4 ++-- .../boundary/open_boundary/method_of_characteristics.jl | 7 +++++-- src/schemes/boundary/open_boundary/system.jl | 9 ++++----- src/visualization/write2vtk.jl | 1 - 4 files changed, 11 insertions(+), 10 deletions(-) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index a52f37673..af558133f 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/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, diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 39613f4b1..b581c0047 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/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) + # 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) # Loop over all fluid neighbors within the kernel cutoff foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 5704b66cf..855d7692e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/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 @@ -367,3 +365,4 @@ end @inline viscosity_model(system::OpenBoundarySPHSystem, neighbor_system::BoundarySystem) = neighbor_system.boundary_model.viscosity # When the neighbor is an open boundary system, just use the viscosity of the fluid `system` instead @inline viscosity_model(system, neighbor_system::OpenBoundarySPHSystem) = system.viscosity +@inline system_sound_speed(system::OpenBoundarySPHSystem) = system_sound_speed(system.fluid_system) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index dabd463dc..e933abfe2 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -325,7 +325,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) From a01311c4ba1fea669dc554613fc17ec6742541c1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Aug 2024 17:34:54 +0200 Subject: [PATCH 06/34] skip empty system --- src/visualization/write2vtk.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index dabd463dc..caac58528 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -81,6 +81,11 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre custom_quantities...) mkpath(output_directory) + # Skip empty systems + if nparticles(system_) == 0 + return + end + # Transfer to CPU if data is on the GPU. Do nothing if already on CPU. v, u, system = transfer2cpu(v_, u_, system_) From 24a21fc1c72923b50d634d8f939259c389290d9d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 Aug 2024 14:29:46 +0200 Subject: [PATCH 07/34] fix test --- test/general/buffer.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/general/buffer.jl b/test/general/buffer.jl index a4117bd61..841945784 100644 --- a/test/general/buffer.jl +++ b/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(), From dcecbe05e37af24531ad242aebd4483c6d58155d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 Aug 2024 14:31:58 +0200 Subject: [PATCH 08/34] fix tests --- test/systems/open_boundary_system.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index 9223a1ef7..beb54b320 100644 --- a/test/systems/open_boundary_system.jl +++ b/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, From 6151f7f243b9f9c2b34231a8ffec857bfcefbeff Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 Aug 2024 14:33:53 +0200 Subject: [PATCH 09/34] fix tests --- test/schemes/boundary/open_boundary/characteristic_variables.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 9636dbfc0..3760f73a8 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/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, From 43e89fd26fb1e4af6e3f4044f620757b92323bec Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 22 Aug 2024 00:09:26 +0200 Subject: [PATCH 10/34] fix bug --- .../density_diffusion.jl | 30 +++++++------------ .../fluid/weakly_compressible_sph/rhs.jl | 11 +++++-- 2 files changed, 19 insertions(+), 22 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index e20951642..2cd92b05a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -51,8 +51,7 @@ struct DensityDiffusionMolteniColagrossi{ELTYPE} <: DensityDiffusion end @inline function density_diffusion_psi(::DensityDiffusionMolteniColagrossi, rho_a, rho_b, - pos_diff, distance, system, neighbor_system, - particle, neighbor) + pos_diff, distance, system, particle, neighbor) return 2 * (rho_a - rho_b) * pos_diff / distance^2 end @@ -87,8 +86,7 @@ struct DensityDiffusionFerrari <: DensityDiffusion end @inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b, - pos_diff, distance, system, neighbor_system, - particle, neighbor) + pos_diff, distance, system, particle, neighbor) (; smoothing_length) = system return ((rho_a - rho_b) / 2smoothing_length) * pos_diff / distance @@ -173,12 +171,11 @@ end @inline function density_diffusion_psi(density_diffusion::DensityDiffusionAntuono, rho_a, rho_b, pos_diff, distance, system, - neighbor_system, particle, neighbor) + particle, neighbor) (; normalized_density_gradient) = density_diffusion normalized_gradient_a = extract_svector(normalized_density_gradient, system, particle) - normalized_gradient_b = extract_svector(normalized_density_gradient, neighbor_system, - neighbor) + normalized_gradient_b = extract_svector(normalized_density_gradient, system, neighbor) # Fist term by Molteni & Colagrossi result = 2 * (rho_a - rho_b) @@ -229,12 +226,9 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search end @inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::FluidSystem, - neighbor_system::FluidSystem, - grad_kernel) + v_particle_system, particle, neighbor, pos_diff, + distance, m_b, rho_a, rho_b, + particle_system::FluidSystem, grad_kernel) # Density diffusion terms are all zero for distance zero distance < sqrt(eps()) && return @@ -245,17 +239,15 @@ end volume_b = m_b / rho_b psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance, - particle_system, neighbor_system, particle, neighbor) + particle_system, particle, neighbor) density_diffusion_term = dot(psi, grad_kernel) * volume_b dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term end # Density diffusion `nothing` or interaction other than fluid-fluid -@inline function density_diffusion!(dv, density_diffusion, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) +@inline function density_diffusion!(dv, density_diffusion, v_particle_system, particle, + neighbor, pos_diff, distance, m_b, rho_a, rho_b, + particle_system, grad_kernel) return dv end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index a1a4f9a37..9e40eacfa 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -111,9 +111,14 @@ end dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel) - density_diffusion!(dv, density_diffusion, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) + # Artificial density diffusion should only be applied to system(s) representing a fluid + # with the same physical properties i.e. density and viscosity. + # TODO: shouldn't be applied to particles on the interface (depends on PR #539) + if particle_system === neighbor_system + density_diffusion!(dv, density_diffusion, v_particle_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, particle_system, + grad_kernel) + end end @inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system, From 0091cc317c755556f6dcfa24383ab6aadf7a1c90 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 22 Aug 2024 12:30:27 +0200 Subject: [PATCH 11/34] fix --- src/general/buffer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/general/buffer.jl b/src/general/buffer.jl index b12cf05f3..cb21a8695 100644 --- a/src/general/buffer.jl +++ b/src/general/buffer.jl @@ -4,7 +4,7 @@ struct SystemBuffer{V} buffer_size :: Int function SystemBuffer(active_size, buffer_size::Integer) - active_particle = vcat(trues(active_size), falses(buffer_size)) + active_particle = vcat(trues(active_size), fill(false, buffer_size)) eachparticle = collect(1:active_size) return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size) From fbaf1082b4108f2a1799bafdadc20ce9e78f749d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 23 Aug 2024 11:27:30 +0200 Subject: [PATCH 12/34] check dimensionality of reference functions --- examples/fluid/pipe_flow_3d.jl | 139 +++++++++++++++++++ src/schemes/boundary/open_boundary/system.jl | 12 ++ 2 files changed, 151 insertions(+) create mode 100644 examples/fluid/pipe_flow_3d.jl diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl new file mode 100644 index 000000000..19a830ce4 --- /dev/null +++ b/examples/fluid/pipe_flow_3d.jl @@ -0,0 +1,139 @@ +# 2D channel flow simulation with open boundaries. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== 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) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.4, 0.4) + +flow_direction = [1.0, 0.0, 0.0] +reynolds_number = 100 +const prescribed_velocity = 2.0 + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) + +fluid_density = 1000.0 + +# For this particular example, it is necessary to have a background pressure. +# Otherwise the suction at the outflow is to big and the simulation becomes unstable. +pressure = 1000.0 + +sound_speed = 10 * prescribed_velocity + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=pressure) + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true, true, true)) + +# 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]^2 + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + +# Alternatively the WCSPH scheme can be used +# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== 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, 0.0) +end + +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) + +open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, + boundary_model=BoundaryModelLastiwka(), + buffer_size=n_buffer_particles, + reference_density=fluid_density, + reference_pressure=pressure, + reference_velocity=velocity_function) + +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) + +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=BoundaryModelLastiwka(), + buffer_size=n_buffer_particles, + reference_density=fluid_density, + reference_pressure=pressure, + reference_velocity=velocity_function) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, + boundary_system) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 855d7692e..b9c4daccb 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -77,6 +77,10 @@ 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 + test_result = reference_velocity(zeros(NDIMS), 0.0) + if length(test_result) != NDIMS + throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) + end reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) end @@ -86,6 +90,10 @@ 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 + test_result = reference_pressure(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_pressure` function must be a scalar function")) + end reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) end @@ -95,6 +103,10 @@ 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 + test_result = reference_density(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_density` function must be a scalar function")) + end reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end From e720529b088a9556ac8e42d03891bf129740a09c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 23 Aug 2024 13:41:09 +0200 Subject: [PATCH 13/34] propagate characteristics --- .../method_of_characteristics.jl | 27 ++++++++++++++++-- src/schemes/boundary/open_boundary/system.jl | 28 +++++++++++-------- 2 files changed, 41 insertions(+), 14 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index b581c0047..152d650e1 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -8,7 +8,8 @@ about the method see [description below](@ref method_of_characteristics). """ struct BoundaryModelLastiwka end -@inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka, +# Called from update callback +@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,19 +46,39 @@ 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 end + +function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) + (; cache) = system + (; previous_characteristics_inited) = cache + + # Propagate characteristics through the open boundary + if !previous_characteristics_inited[1] + + for iteration in 1:4 + calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) + end + + previous_characteristics_inited[1] = true + else + calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) + end + +end + # ==== Characteristics # J1: Associated with convection and entropy and propagates at flow velocity. # J2: Propagates downstream to the local flow # J3: Propagates upstream to the local flow -function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) +function calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) (; volume, cache, boundary_zone) = system - (; characteristics, previous_characteristics) = cache + (; characteristics, previous_characteristics, previous_characteristics_inited) = cache for particle in eachparticle(system) previous_characteristics[1, particle] = characteristics[1, particle] diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index b9c4daccb..6acc41133 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -77,9 +77,11 @@ 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 - test_result = reference_velocity(zeros(NDIMS), 0.0) - if length(test_result) != NDIMS - throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) + 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 @@ -90,9 +92,11 @@ 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 - test_result = reference_pressure(zeros(NDIMS), 0.0) - if length(test_result) != 1 - throw(ArgumentError("`reference_pressure` function must be a scalar function")) + 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 @@ -103,9 +107,11 @@ 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 - test_result = reference_density(zeros(NDIMS), 0.0) - if length(test_result) != 1 - throw(ArgumentError("`reference_density` function must be a scalar function")) + 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 @@ -132,7 +138,7 @@ function create_cache_open_boundary(boundary_model, initial_condition) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) return (; characteristics=characteristics, - previous_characteristics=previous_characteristics) + previous_characteristics=previous_characteristics, previous_characteristics_inited=[false]) end timer_name(::OpenBoundarySPHSystem) = "open_boundary" @@ -202,7 +208,7 @@ 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, + @trixi_timeit timer() "update quantities" update_boundary_quantities!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) From 74004876f044af25e8bf1872aeafbbc106c59201 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 27 Sep 2024 11:55:53 +0200 Subject: [PATCH 14/34] update --- benchmark_results.csv | 7 + benchmark_results_1.csv | 4 + examples/fluid/emptying_tank_2d.jl | 168 ++++++++++++++++++ examples/fluid/filling_tank_2d.jl | 166 +++++++++++++++++ .../method_of_characteristics.jl | 5 +- src/schemes/boundary/open_boundary/system.jl | 10 +- 6 files changed, 352 insertions(+), 8 deletions(-) create mode 100644 benchmark_results.csv create mode 100644 benchmark_results_1.csv create mode 100644 examples/fluid/emptying_tank_2d.jl create mode 100644 examples/fluid/filling_tank_2d.jl diff --git a/benchmark_results.csv b/benchmark_results.csv new file mode 100644 index 000000000..a291a1632 --- /dev/null +++ b/benchmark_results.csv @@ -0,0 +1,7 @@ +Threads,Time +[1],[s] +1,289.78013012 +2,167.728885369 +4,101.710236374 +8,62.295408277 +16,42.730205497 diff --git a/benchmark_results_1.csv b/benchmark_results_1.csv new file mode 100644 index 000000000..b16b75b30 --- /dev/null +++ b/benchmark_results_1.csv @@ -0,0 +1,4 @@ +Threads,Time_s +4,26.705110892 +8,15.746307686 +16,12.644733374 diff --git a/examples/fluid/emptying_tank_2d.jl b/examples/fluid/emptying_tank_2d.jl new file mode 100644 index 000000000..de99f84dd --- /dev/null +++ b/examples/fluid/emptying_tank_2d.jl @@ -0,0 +1,168 @@ +# 2D channel flow simulation with open boundaries. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# 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 +# TODO: filling Tank is not working since the velocity is set. +tspan = (0.0, 4.0) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.4) + +flow_direction = [-1.0, 0.0] +reynolds_number = 100 +const prescribed_velocity = 1.0 + +boundary_size = (domain_size[1], + domain_size[2]) + +fluid_density = 1000.0 + +# For this particular example, it is necessary to have a background pressure. +# Otherwise the suction at the outflow is to big and the simulation becomes unstable. +pressure = 1000.0 + +sound_speed = 10 * prescribed_velocity + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=pressure) + +pipe = RectangularTank(particle_spacing, domain_size.+((boundary_layers+1)*particle_spacing,0.0), boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true)) + +tank = RectangularTank(particle_spacing, (1.5-boundary_layers*particle_spacing, 1.0), (1.5, 1.0), fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, true, true, true), + min_coordinates=(1.0,-0.6)) + +tank_wall_left = RectangularShape(particle_spacing, + (boundary_layers, + round(Int, (1.0 - domain_size[2]) / particle_spacing)), + (1.0, -0.6), density=fluid_density) + + +tank.fluid.coordinates[1, :] .+= particle_spacing * boundary_layers +fluid = union(pipe.fluid, tank.fluid) + +# Shift pipe walls in negative x-direction for the inflow +# pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers + +n_buffer_particles = 200 * pipe.n_particles_per_dimension[2] + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + +# Alternatively the WCSPH scheme can be used +# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +# function velocity_function(pos, t) +# if t > 4 +# return SVector(max((4-t) * prescribed_velocity, 0.0), 0.0) +# else +# return SVector(prescribed_velocity, 0.0) +# end +# 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; fluid_system, +# boundary_model=BoundaryModelLastiwka(), +# buffer_size=n_buffer_particles, +# reference_density=fluid_density, +# reference_pressure=pressure, +# reference_velocity=[prescribed_velocity, 0.0]) + +outflow = OutFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), + flow_direction, open_boundary_layers, density=fluid_density, + particle_spacing) + +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=BoundaryModelLastiwka(), + buffer_size=n_buffer_particles, + reference_density=fluid_density, + reference_pressure=pressure, + reference_velocity=[prescribed_velocity, 0.0]) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +tank_boundary_system = BoundarySPHSystem(tank.boundary, tank_boundary_model) + +wall_boundary_model = BoundaryModelDummyParticles(tank_wall_left.density, + tank_wall_left.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +wall_boundary_system = BoundarySPHSystem(tank_wall_left, wall_boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, open_boundary_out, boundary_system, + tank_boundary_system, wall_boundary_system) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/filling_tank_2d.jl b/examples/fluid/filling_tank_2d.jl new file mode 100644 index 000000000..1654ac401 --- /dev/null +++ b/examples/fluid/filling_tank_2d.jl @@ -0,0 +1,166 @@ +# 2D channel flow simulation with open boundaries. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# 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 +# TODO: filling Tank is not working since the velocity is set. +tspan = (0.0, 4.0) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.4) + +flow_direction = [1.0, 0.0] +reynolds_number = 100 +const prescribed_velocity = 1.0 + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2]) + +fluid_density = 1000.0 + +# For this particular example, it is necessary to have a background pressure. +# Otherwise the suction at the outflow is to big and the simulation becomes unstable. +pressure = 1000.0 + +sound_speed = 10 * prescribed_velocity + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=pressure) + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true)) + +tank = RectangularTank(particle_spacing, (0.0, 0.0), (1.5, 1.0), fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, true, true, true), + min_coordinates=(1.0 + particle_spacing * open_boundary_layers, + -0.6)) + +tank_wall_left = RectangularShape(particle_spacing, + (boundary_layers, + round(Int, (1.0 - domain_size[2]) / particle_spacing)), + (1.0 + particle_spacing * open_boundary_layers, -0.6), + density=fluid_density) + +# Shift pipe walls in negative x-direction for the inflow +pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers + +n_buffer_particles = 200 * pipe.n_particles_per_dimension[2] + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + +# Alternatively the WCSPH scheme can be used +# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +# function velocity_function(pos, t) +# if t > 4 +# return SVector(max((4-t) * prescribed_velocity, 0.0), 0.0) +# else +# return SVector(prescribed_velocity, 0.0) +# end +# 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; fluid_system, + boundary_model=BoundaryModelLastiwka(), + buffer_size=n_buffer_particles, + reference_density=fluid_density, + reference_pressure=pressure, + reference_velocity=[prescribed_velocity, 0.0]) + +# 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; fluid_system, +# boundary_model=BoundaryModelLastiwka(), +# buffer_size=n_buffer_particles, +# reference_density=fluid_density, +# reference_pressure=pressure, +# reference_velocity=velocity_function) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +tank_boundary_system = BoundarySPHSystem(tank.boundary, tank_boundary_model) + +wall_boundary_model = BoundaryModelDummyParticles(tank_wall_left.density, + tank_wall_left.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +wall_boundary_system = BoundarySPHSystem(tank_wall_left, wall_boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, open_boundary_in, boundary_system, + tank_boundary_system, wall_boundary_system) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 152d650e1..3f9d99d0c 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -10,7 +10,7 @@ struct BoundaryModelLastiwka end # Called from update callback @inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, - v, u, v_ode, u_ode, semi, t) + v, u, v_ode, u_ode, semi, t) (; density, pressure, cache, flow_direction, reference_velocity, reference_pressure, reference_density) = system @@ -53,14 +53,12 @@ function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi end end - function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) (; cache) = system (; previous_characteristics_inited) = cache # Propagate characteristics through the open boundary if !previous_characteristics_inited[1] - for iteration in 1:4 calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) end @@ -69,7 +67,6 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) else calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) end - end # ==== Characteristics diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 6acc41133..633ba31b5 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -138,7 +138,8 @@ function create_cache_open_boundary(boundary_model, initial_condition) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) return (; characteristics=characteristics, - previous_characteristics=previous_characteristics, previous_characteristics_inited=[false]) + previous_characteristics=previous_characteristics, + previous_characteristics_inited=[false]) end timer_name(::OpenBoundarySPHSystem) = "open_boundary" @@ -209,9 +210,10 @@ 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_boundary_quantities!(system, - system.boundary_model, - v, u, v_ode, u_ode, - semi, t) + 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) From 550146fd6133c2d3a63fbf2bd07982aab0767a21 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 8 Oct 2024 13:29:57 +0200 Subject: [PATCH 15/34] cleanup --- benchmark_results.csv | 7 ------- benchmark_results_1.csv | 4 ---- 2 files changed, 11 deletions(-) delete mode 100644 benchmark_results.csv delete mode 100644 benchmark_results_1.csv diff --git a/benchmark_results.csv b/benchmark_results.csv deleted file mode 100644 index a291a1632..000000000 --- a/benchmark_results.csv +++ /dev/null @@ -1,7 +0,0 @@ -Threads,Time -[1],[s] -1,289.78013012 -2,167.728885369 -4,101.710236374 -8,62.295408277 -16,42.730205497 diff --git a/benchmark_results_1.csv b/benchmark_results_1.csv deleted file mode 100644 index b16b75b30..000000000 --- a/benchmark_results_1.csv +++ /dev/null @@ -1,4 +0,0 @@ -Threads,Time_s -4,26.705110892 -8,15.746307686 -16,12.644733374 From 90766fea9ca0a07ed1a100d2bf29016d14f53ed4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 01:33:33 +0200 Subject: [PATCH 16/34] update --- .../boundary/open_boundary/method_of_characteristics.jl | 2 +- src/schemes/boundary/open_boundary/system.jl | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 44f8ce7c1..d23d66fa9 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -8,7 +8,7 @@ about the method see [description below](@ref method_of_characteristics). """ struct BoundaryModelLastiwka end -# Called from update callback +# 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, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 633ba31b5..f846c1c6e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -209,11 +209,10 @@ 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_boundary_quantities!(system, + @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, system.boundary_model, v, u, v_ode, - u_ode, - semi, t) + u_ode, semi, t) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) @@ -385,4 +384,3 @@ end @inline viscosity_model(system::OpenBoundarySPHSystem, neighbor_system::BoundarySystem) = neighbor_system.boundary_model.viscosity # When the neighbor is an open boundary system, just use the viscosity of the fluid `system` instead @inline viscosity_model(system, neighbor_system::OpenBoundarySPHSystem) = system.viscosity -@inline system_sound_speed(system::OpenBoundarySPHSystem) = system_sound_speed(system.fluid_system) From 4ed0b91a1621d8e022a13f1bb2a62c78aa80c483 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 13:02:07 +0200 Subject: [PATCH 17/34] Increase errors for 1.11 --- test/validation/validation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/validation/validation.jl b/test/validation/validation.jl index 1ee36f1bd..495d14de7 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -47,13 +47,14 @@ @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 - if VERSION >= v"1.10" + if VERSION == v"1.10" @test isapprox(error_edac_P1, 0, atol=eps()) @test isapprox(error_edac_P2, 0, atol=eps()) @test isapprox(error_wcsph_P1, 0, atol=eps()) @test isapprox(error_wcsph_P2, 0, atol=eps()) else # 1.9 causes a large difference in the solution + # 1.11 requires a performance hotfix which will likely change these results again @test isapprox(error_edac_P1, 0, atol=4e-9) @test isapprox(error_edac_P2, 0, atol=3e-11) @test isapprox(error_wcsph_P1, 0, atol=26.3) From 1ce04600e8aa846e69e130bca2b1be84d7c95eff Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 13:09:12 +0200 Subject: [PATCH 18/34] Fix invalidations --- .github/workflows/Invalidations.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 9b2831d8c..e2116321e 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -18,7 +18,8 @@ jobs: steps: - uses: julia-actions/setup-julia@v2 with: - version: '1' + # TODO: Is broken in 1.11 revert to 1 after fix + version: '1.10' - uses: actions/checkout@v4 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-invalidations@v1 From 1774f5a1b81ea3ac5e4ab60c94ef0572cb35ad35 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 14:52:10 +0200 Subject: [PATCH 19/34] Fix tests --- .github/workflows/ci.yml | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 04070512f..252cae3a3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -27,7 +27,6 @@ on: - 'docs/**' workflow_dispatch: - # Cancel redundant CI tests automatically concurrency: group: ${{ github.workflow }}-${{ github.ref }} @@ -36,25 +35,26 @@ concurrency: jobs: build: name: Run Tests (Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}) - runs-on: ubuntu-latest + runs-on: ${{ matrix.os }} strategy: - # Don't cancel all running jobs when one job fails fail-fast: false matrix: - version: - - '1.9' - - '1.10' - - '1' - os: - - ubuntu-latest - arch: - - x64 include: - # Also run tests on Windows and macOS-ARM, but only with the latest Julia version - # TODO: change back to '1' after JuliaLang/julia#55009 has been backported + # Ubuntu jobs with Julia 1.9, 1.10, and 1.11 + - version: '1.9' + os: ubuntu-latest + arch: x64 + - version: '1.10' + os: ubuntu-latest + arch: x64 + - version: '1' + os: ubuntu-latest + arch: x64 + # Windows job with Julia 1.10 - version: '1.10' os: windows-latest arch: x64 + # macOS ARM job with Julia 1.10 - version: '1.10' os: macos-14 arch: arm64 @@ -62,29 +62,35 @@ jobs: steps: - name: Check out project uses: actions/checkout@v4 + - name: Set up Julia uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} + - uses: julia-actions/cache@v2 + - name: Build package uses: julia-actions/julia-buildpkg@v1 + - name: Run unit tests uses: julia-actions/julia-runtest@v1 with: annotate: true - # Only run coverage in one Job (Ubuntu and latest Julia version) + # Only run coverage in the Ubuntu job with Julia 1.11 coverage: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' }} env: TRIXIPARTICLES_TEST: unit + - name: Process coverage results - # Only run coverage in one Job (Ubuntu and latest Julia version) + # Only process coverage in the Ubuntu job with Julia 1.11 if: matrix.os == 'ubuntu-latest' && matrix.version == '1' uses: julia-actions/julia-processcoverage@v1 with: directories: src,test + - name: Upload coverage report to Codecov - # Only run coverage in one Job (Ubuntu and latest Julia version) + # Only upload coverage in the Ubuntu job with Julia 1.11 if: matrix.os == 'ubuntu-latest' && matrix.version == '1' uses: codecov/codecov-action@v4 with: @@ -93,6 +99,7 @@ jobs: flags: unit env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + - name: Run example tests uses: julia-actions/julia-runtest@v1 with: From 4daf98407729ee0e79af8973136d0b307cb8cf52 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 14:54:02 +0200 Subject: [PATCH 20/34] Update ci.yml --- .github/workflows/ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 252cae3a3..2241ca9e4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -77,20 +77,20 @@ jobs: uses: julia-actions/julia-runtest@v1 with: annotate: true - # Only run coverage in the Ubuntu job with Julia 1.11 + # Only run coverage in one Job (Ubuntu and latest Julia version) coverage: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' }} env: TRIXIPARTICLES_TEST: unit - name: Process coverage results - # Only process coverage in the Ubuntu job with Julia 1.11 + # Only run coverage in one Job (Ubuntu and latest Julia version) if: matrix.os == 'ubuntu-latest' && matrix.version == '1' uses: julia-actions/julia-processcoverage@v1 with: directories: src,test - name: Upload coverage report to Codecov - # Only upload coverage in the Ubuntu job with Julia 1.11 + # Only run coverage in one Job (Ubuntu and latest Julia version) if: matrix.os == 'ubuntu-latest' && matrix.version == '1' uses: codecov/codecov-action@v4 with: From ac2eb2c93a28701affffdec5bf18f7edd2bdd62f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 15:21:26 +0200 Subject: [PATCH 21/34] revert --- .github/workflows/ci.yml | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2241ca9e4..6091cc84d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -35,29 +35,28 @@ concurrency: jobs: build: name: Run Tests (Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}) - runs-on: ${{ matrix.os }} + runs-on: ubuntu-latest strategy: fail-fast: false matrix: - include: - # Ubuntu jobs with Julia 1.9, 1.10, and 1.11 - - version: '1.9' - os: ubuntu-latest - arch: x64 - - version: '1.10' - os: ubuntu-latest - arch: x64 - - version: '1' - os: ubuntu-latest - arch: x64 - # Windows job with Julia 1.10 - - version: '1.10' - os: windows-latest - arch: x64 - # macOS ARM job with Julia 1.10 - - version: '1.10' - os: macos-14 - arch: arm64 + version: + - '1.9' + - '1.10' + - '1' + os: + - ubuntu-latest + arch: + - x64 + include: + # Also run tests on Windows and macOS-ARM, but only with the latest Julia version + # TODO: change back to '1' after JuliaLang/julia#55009 has been backported + - version: '1.10' + os: windows-latest + arch: x64 + # macOS ARM job with Julia 1.10 + - version: '1.10' + os: macos-14 + arch: arm64 steps: - name: Check out project From 258f2baadbb6caca916f039f53096da4a7626e8f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Oct 2024 15:28:35 +0200 Subject: [PATCH 22/34] Update ci.yml --- .github/workflows/ci.yml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6091cc84d..06c868258 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -48,15 +48,19 @@ jobs: arch: - x64 include: - # Also run tests on Windows and macOS-ARM, but only with the latest Julia version - # TODO: change back to '1' after JuliaLang/julia#55009 has been backported + # Also run tests on Windows and macOS-ARM, but only with the latest Julia version and 1.10 - version: '1.10' os: windows-latest arch: x64 - # macOS ARM job with Julia 1.10 - version: '1.10' os: macos-14 arch: arm64 + - version: '1' + os: windows-latest + arch: x64 + - version: '1' + os: macos-14 + arch: arm64 steps: - name: Check out project From 1a9f3f5b1ec5422a9445e906d8f7fe4dc353a4b1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 11:41:01 +0200 Subject: [PATCH 23/34] Update test/validation/validation.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- test/validation/validation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/validation/validation.jl b/test/validation/validation.jl index 495d14de7..01efa6b6a 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -54,7 +54,7 @@ @test isapprox(error_wcsph_P2, 0, atol=eps()) else # 1.9 causes a large difference in the solution - # 1.11 requires a performance hotfix which will likely change these results again + # TODO 1.11 requires a performance hotfix which will likely change these results again @test isapprox(error_edac_P1, 0, atol=4e-9) @test isapprox(error_edac_P2, 0, atol=3e-11) @test isapprox(error_wcsph_P1, 0, atol=26.3) From da9be38a4c23f747d8fb43e3f0b0b373cb229622 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 12:41:00 +0200 Subject: [PATCH 24/34] revert changes that had no benefit --- .../method_of_characteristics.jl | 20 ++----------------- src/schemes/boundary/open_boundary/system.jl | 3 +-- 2 files changed, 3 insertions(+), 20 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index d23d66fa9..344222f13 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -53,29 +53,13 @@ function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi end end -function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) - (; cache) = system - (; previous_characteristics_inited) = cache - - # Propagate characteristics through the open boundary - if !previous_characteristics_inited[1] - for iteration in 1:4 - calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) - end - - previous_characteristics_inited[1] = true - else - calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) - end -end - # ==== Characteristics # J1: Associated with convection and entropy and propagates at flow velocity. # J2: Propagates downstream to the local flow # J3: Propagates upstream to the local flow -function calc_characteristics!(system, v, u, v_ode, u_ode, semi, t) +function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) (; volume, cache, boundary_zone) = system - (; characteristics, previous_characteristics, previous_characteristics_inited) = cache + (; characteristics, previous_characteristics) = cache for particle in eachparticle(system) previous_characteristics[1, particle] = characteristics[1, particle] diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index f846c1c6e..c62ac6475 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -138,8 +138,7 @@ function create_cache_open_boundary(boundary_model, initial_condition) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) return (; characteristics=characteristics, - previous_characteristics=previous_characteristics, - previous_characteristics_inited=[false]) + previous_characteristics=previous_characteristics) end timer_name(::OpenBoundarySPHSystem) = "open_boundary" From 95be0e9077a30dacabc387e37993f3f93c0f5a36 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 12:52:58 +0200 Subject: [PATCH 25/34] update --- .../boundary/open_boundary/boundary_zones.jl | 5 +++-- .../open_boundary/method_of_characteristics.jl | 17 +++++++++++++---- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index cc9a77157..098436ddd 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -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(dot(normalize(edge1), normalize(edge2)), 1.0; atol=1e-7) + 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 344222f13..1a5e0b083 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -51,6 +51,8 @@ function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end + + return system end # ==== Characteristics @@ -100,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 avg_J1 + avg_J2 + avg_J3 > eps() + 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] @@ -166,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 From 4b87c7c844a85228594e26848245b8c615e50d92 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 13:48:10 +0200 Subject: [PATCH 26/34] cleanup --- examples/fluid/emptying_tank_2d.jl | 168 ------------------- examples/fluid/filling_tank_2d.jl | 166 ------------------ examples/fluid/pipe_flow_2d.jl | 6 +- examples/fluid/pipe_flow_3d.jl | 155 ++++------------- src/schemes/boundary/open_boundary/system.jl | 8 +- 5 files changed, 38 insertions(+), 465 deletions(-) delete mode 100644 examples/fluid/emptying_tank_2d.jl delete mode 100644 examples/fluid/filling_tank_2d.jl diff --git a/examples/fluid/emptying_tank_2d.jl b/examples/fluid/emptying_tank_2d.jl deleted file mode 100644 index de99f84dd..000000000 --- a/examples/fluid/emptying_tank_2d.jl +++ /dev/null @@ -1,168 +0,0 @@ -# 2D channel flow simulation with open boundaries. - -using TrixiParticles -using OrdinaryDiffEq - -# ========================================================================================== -# ==== Resolution -particle_spacing = 0.05 - -# Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 4 - -# 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 -# TODO: filling Tank is not working since the velocity is set. -tspan = (0.0, 4.0) - -# Boundary geometry and initial fluid particle positions -domain_size = (1.0, 0.4) - -flow_direction = [-1.0, 0.0] -reynolds_number = 100 -const prescribed_velocity = 1.0 - -boundary_size = (domain_size[1], - domain_size[2]) - -fluid_density = 1000.0 - -# For this particular example, it is necessary to have a background pressure. -# Otherwise the suction at the outflow is to big and the simulation becomes unstable. -pressure = 1000.0 - -sound_speed = 10 * prescribed_velocity - -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, background_pressure=pressure) - -pipe = RectangularTank(particle_spacing, domain_size.+((boundary_layers+1)*particle_spacing,0.0), boundary_size, fluid_density, - pressure=pressure, n_layers=boundary_layers, - faces=(false, false, true, true)) - -tank = RectangularTank(particle_spacing, (1.5-boundary_layers*particle_spacing, 1.0), (1.5, 1.0), fluid_density, - pressure=pressure, n_layers=boundary_layers, - faces=(false, true, true, true), - min_coordinates=(1.0,-0.6)) - -tank_wall_left = RectangularShape(particle_spacing, - (boundary_layers, - round(Int, (1.0 - domain_size[2]) / particle_spacing)), - (1.0, -0.6), density=fluid_density) - - -tank.fluid.coordinates[1, :] .+= particle_spacing * boundary_layers -fluid = union(pipe.fluid, tank.fluid) - -# Shift pipe walls in negative x-direction for the inflow -# pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers - -n_buffer_particles = 200 * pipe.n_particles_per_dimension[2] - -# ========================================================================================== -# ==== Fluid -smoothing_length = 3.0 * particle_spacing -smoothing_kernel = WendlandC2Kernel{2}() - -fluid_density_calculator = ContinuityDensity() - -kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number - -viscosity = ViscosityAdami(nu=kinematic_viscosity) - -fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=viscosity, - density_calculator=fluid_density_calculator, - buffer_size=n_buffer_particles) - -# Alternatively the WCSPH scheme can be used -# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) -# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) - -# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, -# state_equation, smoothing_kernel, -# smoothing_length, viscosity=viscosity, -# buffer_size=n_buffer_particles) - -# ========================================================================================== -# ==== Open Boundary -# function velocity_function(pos, t) -# if t > 4 -# return SVector(max((4-t) * prescribed_velocity, 0.0), 0.0) -# else -# return SVector(prescribed_velocity, 0.0) -# end -# 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; fluid_system, -# boundary_model=BoundaryModelLastiwka(), -# buffer_size=n_buffer_particles, -# reference_density=fluid_density, -# reference_pressure=pressure, -# reference_velocity=[prescribed_velocity, 0.0]) - -outflow = OutFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), - flow_direction, open_boundary_layers, density=fluid_density, - particle_spacing) - -open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=BoundaryModelLastiwka(), - buffer_size=n_buffer_particles, - reference_density=fluid_density, - reference_pressure=pressure, - reference_velocity=[prescribed_velocity, 0.0]) - -# ========================================================================================== -# ==== Boundary - -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) - -tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -tank_boundary_system = BoundarySPHSystem(tank.boundary, tank_boundary_model) - -wall_boundary_model = BoundaryModelDummyParticles(tank_wall_left.density, - tank_wall_left.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -wall_boundary_system = BoundarySPHSystem(tank_wall_left, wall_boundary_model) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, open_boundary_out, boundary_system, - tank_boundary_system, wall_boundary_system) - -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") - -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) - -sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing - save_everystep=false, callback=callbacks); diff --git a/examples/fluid/filling_tank_2d.jl b/examples/fluid/filling_tank_2d.jl deleted file mode 100644 index 1654ac401..000000000 --- a/examples/fluid/filling_tank_2d.jl +++ /dev/null @@ -1,166 +0,0 @@ -# 2D channel flow simulation with open boundaries. - -using TrixiParticles -using OrdinaryDiffEq - -# ========================================================================================== -# ==== Resolution -particle_spacing = 0.05 - -# Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 4 - -# 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 -# TODO: filling Tank is not working since the velocity is set. -tspan = (0.0, 4.0) - -# Boundary geometry and initial fluid particle positions -domain_size = (1.0, 0.4) - -flow_direction = [1.0, 0.0] -reynolds_number = 100 -const prescribed_velocity = 1.0 - -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2]) - -fluid_density = 1000.0 - -# For this particular example, it is necessary to have a background pressure. -# Otherwise the suction at the outflow is to big and the simulation becomes unstable. -pressure = 1000.0 - -sound_speed = 10 * prescribed_velocity - -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, background_pressure=pressure) - -pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, - pressure=pressure, n_layers=boundary_layers, - faces=(false, false, true, true)) - -tank = RectangularTank(particle_spacing, (0.0, 0.0), (1.5, 1.0), fluid_density, - pressure=pressure, n_layers=boundary_layers, - faces=(false, true, true, true), - min_coordinates=(1.0 + particle_spacing * open_boundary_layers, - -0.6)) - -tank_wall_left = RectangularShape(particle_spacing, - (boundary_layers, - round(Int, (1.0 - domain_size[2]) / particle_spacing)), - (1.0 + particle_spacing * open_boundary_layers, -0.6), - density=fluid_density) - -# Shift pipe walls in negative x-direction for the inflow -pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers - -n_buffer_particles = 200 * pipe.n_particles_per_dimension[2] - -# ========================================================================================== -# ==== Fluid -smoothing_length = 3.0 * particle_spacing -smoothing_kernel = WendlandC2Kernel{2}() - -fluid_density_calculator = ContinuityDensity() - -kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number - -viscosity = ViscosityAdami(nu=kinematic_viscosity) - -fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=viscosity, - density_calculator=fluid_density_calculator, - buffer_size=n_buffer_particles) - -# Alternatively the WCSPH scheme can be used -# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) -# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) - -# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, -# state_equation, smoothing_kernel, -# smoothing_length, viscosity=viscosity, -# buffer_size=n_buffer_particles) - -# ========================================================================================== -# ==== Open Boundary -# function velocity_function(pos, t) -# if t > 4 -# return SVector(max((4-t) * prescribed_velocity, 0.0), 0.0) -# else -# return SVector(prescribed_velocity, 0.0) -# end -# 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; fluid_system, - boundary_model=BoundaryModelLastiwka(), - buffer_size=n_buffer_particles, - reference_density=fluid_density, - reference_pressure=pressure, - reference_velocity=[prescribed_velocity, 0.0]) - -# 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; fluid_system, -# boundary_model=BoundaryModelLastiwka(), -# buffer_size=n_buffer_particles, -# reference_density=fluid_density, -# reference_pressure=pressure, -# reference_velocity=velocity_function) - -# ========================================================================================== -# ==== Boundary - -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) - -tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -tank_boundary_system = BoundarySPHSystem(tank.boundary, tank_boundary_model) - -wall_boundary_model = BoundaryModelDummyParticles(tank_wall_left.density, - tank_wall_left.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -wall_boundary_system = BoundarySPHSystem(tank_wall_left, wall_boundary_model) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, open_boundary_in, boundary_system, - tank_boundary_system, wall_boundary_system) - -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") - -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) - -sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing - save_everystep=false, callback=callbacks); diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index af558133f..46adefbd0 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/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 diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl index 19a830ce4..5e9994a53 100644 --- a/examples/fluid/pipe_flow_3d.jl +++ b/examples/fluid/pipe_flow_3d.jl @@ -1,139 +1,44 @@ -# 2D channel flow simulation with open boundaries. - +# 3D channel flow simulation with open boundaries. using TrixiParticles using OrdinaryDiffEq -# ========================================================================================== -# ==== 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) - -# Boundary geometry and initial fluid particle positions -domain_size = (1.0, 0.4, 0.4) - -flow_direction = [1.0, 0.0, 0.0] -reynolds_number = 100 -const prescribed_velocity = 2.0 - -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2], domain_size[3]) - -fluid_density = 1000.0 - -# For this particular example, it is necessary to have a background pressure. -# Otherwise the suction at the outflow is to big and the simulation becomes unstable. -pressure = 1000.0 - -sound_speed = 10 * prescribed_velocity - -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, background_pressure=pressure) - -pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, - pressure=pressure, n_layers=boundary_layers, - faces=(false, false, true, true, true, true)) - -# 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]^2 +# load variables +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + sol=nothing) -# ========================================================================================== -# ==== Fluid -smoothing_length = 3.0 * particle_spacing -smoothing_kernel = WendlandC2Kernel{3}() - -fluid_density_calculator = ContinuityDensity() - -kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number - -viscosity = ViscosityAdami(nu=kinematic_viscosity) - -fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=viscosity, - density_calculator=fluid_density_calculator, - buffer_size=n_buffer_particles) - -# Alternatively the WCSPH scheme can be used -# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) -# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) - -# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, -# state_equation, smoothing_kernel, -# smoothing_length, viscosity=viscosity, -# buffer_size=n_buffer_particles) - -# ========================================================================================== -# ==== Open Boundary -function velocity_function(pos, t) +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 -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) - -open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelLastiwka(), - buffer_size=n_buffer_particles, - reference_density=fluid_density, - reference_pressure=pressure, - reference_velocity=velocity_function) - -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) - -open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=BoundaryModelLastiwka(), - buffer_size=n_buffer_particles, - reference_density=fluid_density, - reference_pressure=pressure, - reference_velocity=velocity_function) - -# ========================================================================================== -# ==== Boundary - -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), - smoothing_kernel, smoothing_length) - -boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, - boundary_system) +domain_size = (1.0, 0.4, 0.4) -ode = semidiscretize(semi, tspan) +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) -info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") +pipe3d = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true, true, true)) -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) +flow_direction = [1.0, 0.0, 0.0] -sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing - save_everystep=false, callback=callbacks); +# 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)) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index c62ac6475..bad41673f 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -209,9 +209,11 @@ 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 boundary quantities" update_boundary_quantities!(system, - system.boundary_model, - v, u, v_ode, - u_ode, semi, t) + 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) From 8c7257643859c9d295ce432d5533cdf1d698a7f1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 13:48:49 +0200 Subject: [PATCH 27/34] include in test run --- test/examples/examples.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index c962b5912..81d5a608a 100644 --- a/test/examples/examples.jl +++ b/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", From 782f9dedd525122efd226a73eafeb392233c19cf Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 15:32:20 +0200 Subject: [PATCH 28/34] remove redundancy --- .../boundary/open_boundary/boundary_zones.jl | 228 ++++++++---------- 1 file changed, 106 insertions(+), 122 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 098436ddd..f1fffba96 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -1,3 +1,76 @@ +abstract type FlowBoundary{NDIMS} end + +struct InFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD +end + +struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD +end + +function _create_flow_boundary(; plane, flow_direction, density, particle_spacing, + initial_condition, extrude_geometry, + open_boundary_layers::Integer, is_inflow::Bool) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector{length(flow_direction)}(flow_direction)) + + # Determine extrusion direction + extrusion_direction = is_inflow ? -flow_direction_ : flow_direction_ + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=extrusion_direction, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, + density, + direction=extrusion_direction, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_set` is normal to the plane. + # The normal vector must point in the correct direction. + normal_vector = normalize(spanning_set[:, 1]) + dot_ = dot(normal_vector, flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` must be normal to the plane defined by `plane` points.")) + end + + # Flip the normal vector to point in the correct direction + spanning_set[:, 1] .*= is_inflow ? -sign(dot_) : sign(dot_) + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove particles outside the boundary zone + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ +end + @doc raw""" InFlow(; plane, flow_direction, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, @@ -71,67 +144,22 @@ inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, dens !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct InFlow{NDIMS, IC, S, ZO, ZW, FD} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD - - function InFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - if open_boundary_layers <= 0 - throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) - end - - # Unit vector pointing in downstream direction - flow_direction_ = normalize(SVector(flow_direction...)) - - # Sample particles in boundary zone - if isnothing(initial_condition) && isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, - density, - direction=-flow_direction_, - n_extrude=open_boundary_layers) - elseif !isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; - particle_spacing, - density, - direction=-flow_direction_, - n_extrude=open_boundary_layers) - end - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - zone_width = open_boundary_layers * initial_condition.particle_spacing - - # Vectors spanning the boundary zone/box - spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) - - # First vector of `spanning_vectors` is normal to the inflow plane. - # The normal vector must point in upstream direction for an inflow boundary. - dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) - - 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 - - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) - - # Remove particles outside the boundary zone. - # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. - ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) - - return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), - typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, - flow_direction_) - end +function InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(; + plane, + flow_direction, + density, + particle_spacing, + initial_condition, + extrude_geometry, + open_boundary_layers, + is_inflow=true) + return InFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_, + zone_origin, zone_width, + flow_direction_) end @doc raw""" @@ -207,66 +235,22 @@ outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, de !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD - - function OutFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - if open_boundary_layers <= 0 - throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) - end - - # Unit vector pointing in downstream direction - flow_direction_ = normalize(SVector(flow_direction...)) - - # Sample particles in boundary zone - if isnothing(initial_condition) && isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, - density, - direction=flow_direction_, - n_extrude=open_boundary_layers) - elseif !isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; - particle_spacing, density, - direction=flow_direction_, - n_extrude=open_boundary_layers) - end - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - zone_width = open_boundary_layers * initial_condition.particle_spacing - - # Vectors spanning the boundary zone/box - spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) - - # First vector of `spanning_vectors` is normal to the outflow plane. - # The normal vector must point in downstream direction for an outflow boundary. - dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) - - 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 - - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) - - # Remove particles outside the boundary zone. - # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. - ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) - - return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), - typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, - flow_direction_) - end +function OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(; + plane, + flow_direction, + density, + particle_spacing, + initial_condition, + extrude_geometry, + open_boundary_layers, + is_inflow=false) + return OutFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_, + zone_origin, zone_width, + flow_direction_) end @inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS @@ -300,7 +284,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) return hcat(c, edge1, edge2) end -@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords) +@inline function is_in_boundary_zone(boundary_zone::FlowBoundary, particle_coords) (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin From 6f5ca10886b258421c5acecaea1c195737997316 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 10 Oct 2024 16:56:42 +0200 Subject: [PATCH 29/34] revert --- .../boundary/open_boundary/boundary_zones.jl | 229 ++++++++++-------- 1 file changed, 123 insertions(+), 106 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index f1fffba96..03e720166 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -1,76 +1,3 @@ -abstract type FlowBoundary{NDIMS} end - -struct InFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD -end - -struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} <: FlowBoundary{NDIMS} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD -end - -function _create_flow_boundary(; plane, flow_direction, density, particle_spacing, - initial_condition, extrude_geometry, - open_boundary_layers::Integer, is_inflow::Bool) - if open_boundary_layers <= 0 - throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) - end - - # Unit vector pointing in downstream direction - flow_direction_ = normalize(SVector{length(flow_direction)}(flow_direction)) - - # Determine extrusion direction - extrusion_direction = is_inflow ? -flow_direction_ : flow_direction_ - - # Sample particles in boundary zone - if isnothing(initial_condition) && isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, - density, - direction=extrusion_direction, - n_extrude=open_boundary_layers) - elseif !isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; - particle_spacing, - density, - direction=extrusion_direction, - n_extrude=open_boundary_layers) - end - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - zone_width = open_boundary_layers * initial_condition.particle_spacing - - # Vectors spanning the boundary zone/box - spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) - - # First vector of `spanning_set` is normal to the plane. - # The normal vector must point in the correct direction. - normal_vector = normalize(spanning_set[:, 1]) - dot_ = dot(normal_vector, flow_direction_) - - if !isapprox(abs(dot_), 1.0, atol=1e-7) - throw(ArgumentError("`flow_direction` must be normal to the plane defined by `plane` points.")) - end - - # Flip the normal vector to point in the correct direction - spanning_set[:, 1] .*= is_inflow ? -sign(dot_) : sign(dot_) - - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) - - # Remove particles outside the boundary zone - ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) - - return NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ -end - @doc raw""" InFlow(; plane, flow_direction, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, @@ -144,22 +71,67 @@ inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, dens !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -function InFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(; - plane, - flow_direction, - density, - particle_spacing, - initial_condition, - extrude_geometry, - open_boundary_layers, - is_inflow=true) - return InFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_, - zone_origin, zone_width, - flow_direction_) +struct InFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the inflow plane. + # The normal vector must point in upstream direction for an inflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to inflow plane")) + 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. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end end @doc raw""" @@ -235,22 +207,67 @@ outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, de !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -function OutFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - NDIMS, ic, spanning_set_, zone_origin, zone_width, flow_direction_ = _create_flow_boundary(; - plane, - flow_direction, - density, - particle_spacing, - initial_condition, - extrude_geometry, - open_boundary_layers, - is_inflow=false) - return OutFlow{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), typeof(flow_direction_)}(ic, spanning_set_, - zone_origin, zone_width, - flow_direction_) +struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, density, + direction=flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the outflow plane. + # The normal vector must point in downstream direction for an outflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to outflow plane")) + 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. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end end @inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS @@ -284,7 +301,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) return hcat(c, edge1, edge2) end -@inline function is_in_boundary_zone(boundary_zone::FlowBoundary, particle_coords) +@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords) (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin From 5e688d23e6dddc35f3c4d392d0d5cf71cae59467 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 11 Oct 2024 14:49:28 +0200 Subject: [PATCH 30/34] fix tests --- src/schemes/boundary/open_boundary/boundary_zones.jl | 1 - src/schemes/boundary/open_boundary/method_of_characteristics.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 03e720166..24a544562 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -256,7 +256,6 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} # 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. diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 1a5e0b083..9cef5276b 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -103,7 +103,7 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end # To prevent NANs here if the boundary has not been in contact before. - if avg_J1 + avg_J2 + avg_J3 > eps() + if abs(avg_J1 + avg_J2 + avg_J3) > eps() characteristics[1, particle] = avg_J1 / counter characteristics[2, particle] = avg_J2 / counter characteristics[3, particle] = avg_J3 / counter From 7557879e754815e2e6d04348da420fe99cbd27c0 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 11 Oct 2024 14:52:01 +0200 Subject: [PATCH 31/34] fix --- src/schemes/boundary/open_boundary/method_of_characteristics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 9cef5276b..7fe5f0d39 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -103,7 +103,7 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end # To prevent NANs here if the boundary has not been in contact before. - if abs(avg_J1 + avg_J2 + avg_J3) > eps() + if counter > 0 characteristics[1, particle] = avg_J1 / counter characteristics[2, particle] = avg_J2 / counter characteristics[3, particle] = avg_J3 / counter From 8294aa9e587764e28bae5613c5ac2b76494c6a5e Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 11 Oct 2024 17:23:18 +0200 Subject: [PATCH 32/34] fix test --- test/schemes/boundary/open_boundary/boundary_zone.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index 7a34a4bca..afa812fd1 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, From bdd67568e7ab090209d21d4333cdafa98efd8278 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 11 Oct 2024 17:41:38 +0200 Subject: [PATCH 33/34] fix test --- src/schemes/boundary/open_boundary/boundary_zones.jl | 2 +- src/setups/extrude_geometry.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 24a544562..35c2ba212 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -290,7 +290,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) edge2 = plane_points[3] - plane_points[1] # Check if the edges are linearly dependent (to avoid degenerate planes) - if isapprox(dot(normalize(edge1), normalize(edge2)), 1.0; atol=1e-7) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) throw(ArgumentError("The vectors `AB` and `AC` must not be collinear")) end diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index c666ec3d7..0846bbb25 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 From 1ec7eacd399ea1059f3b95b18965e3f874ac01c8 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 16 Oct 2024 13:39:02 +0200 Subject: [PATCH 34/34] fix the test --- examples/fluid/pipe_flow_3d.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl index 5e9994a53..71506c37b 100644 --- a/examples/fluid/pipe_flow_3d.jl +++ b/examples/fluid/pipe_flow_3d.jl @@ -2,6 +2,8 @@ using TrixiParticles using OrdinaryDiffEq +tspan = (0.0, 2.0) + # load variables trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), sol=nothing)