diff --git a/examples/fluid/simple_advection_2d.jl b/examples/fluid/simple_advection_2d.jl new file mode 100644 index 0000000000..5fef28a5ff --- /dev/null +++ b/examples/fluid/simple_advection_2d.jl @@ -0,0 +1,144 @@ +# This is a simple advection example with open boundaries. The setup is similar to +# `periodic_channel_2d.jl`, but uses open boundaries instead of periodic ones. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.02 + +# 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. +open_boundary_layers = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.5) + +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 + +sound_speed = 10 * prescribed_velocity + +state_equation = nothing + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + velocity=[prescribed_velocity, 0.0], + n_layers=boundary_layers, faces=(false, false, true, true)) + +# Shift pipe walls in negative x-direction for the inflow +pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers + +NDIMS = ndims(pipe.fluid) + +n_buffer_particles = 5 * pipe.n_particles_per_dimension[2]^(NDIMS - 1) + +# ========================================================================================== +# ==== Fluid +wcsph = true + +smoothing_length = 1.2 * particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{NDIMS}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +# Alternatively the EDAC scheme can be used +if wcsph + state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) + + fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + buffer_size=n_buffer_particles) +else + fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, + smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) +end + +# ========================================================================================== +# ==== Open Boundary + +velocity_function2d(pos, t) = SVector(prescribed_velocity, 0.0) + +open_boundary_model = BoundaryModelTafuni() + +boundary_type_in = InFlow() +plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers, + density=fluid_density, particle_spacing, + boundary_type=boundary_type_in) + +open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles, + reference_velocity=velocity_function2d) + +boundary_type_out = OutFlow() +plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) +outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + boundary_type=boundary_type_out) +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles, + reference_velocity=velocity_function2d) + +# ========================================================================================== +# ==== Boundary +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2) +max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2) + +nhs = GridNeighborhoodSearch{NDIMS}(; cell_list=FullGridCellList(; min_corner, max_corner), + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, + boundary_system, neighborhood_search=nhs, + parallelization_backend=PolyesterBackend()) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +particle_shifting = nothing + +extra_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), + particle_shifting, extra_callback) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-8, # 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/system.jl b/src/schemes/boundary/open_boundary/system.jl index f299fd9286..e0a2146693 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -91,11 +91,7 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; initial_condition = allocate_buffer(initial_condition, buffer) NDIMS = ndims(initial_condition) - - pressure = copy(initial_condition.pressure) - mass = copy(initial_condition.mass) - density = copy(initial_condition.density) - volume = similar(initial_condition.density) + ELTYPE = eltype(initial_condition) if !(reference_velocity isa Function || isnothing(reference_velocity) || (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) @@ -143,6 +139,21 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end + coordinates_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, + initial_condition.coordinates) + velocities_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, + initial_condition.velocity) + + initial_condition.velocity .= stack(reference_value.(reference_velocity_, + velocities_svector, + coordinates_svector, 0)) + pressure = copy(reference_value.(reference_pressure_, initial_condition.pressure, + coordinates_svector, 0)) + density = copy(reference_value.(reference_density_, initial_condition.density, + coordinates_svector, 0)) + mass = copy(initial_condition.mass) + volume = similar(initial_condition.density) + cache = create_cache_open_boundary(boundary_model, initial_condition, reference_density, reference_velocity, reference_pressure) diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index ede9e6bbf9..eebb273c22 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -289,6 +289,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/simple_advection_2d.jl.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "simple_advection_2d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (WCSPH)" begin @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), joinpath(examples_dir(), "fluid",