Fix IC in open boundaries and add simple_advection_2d.jl

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);

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)

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",