Revise open boundaries (#866)

* swap extrapolation

* swap density pressure

* adapt docs

* rm example

* fix missing nhs update

* implement different mirror methods

* add docs

* fix tests

* fix typos

* add docs

* add tests

* apply formatter

* add `average_velocity!` for `BoundaryModelLastiwka`

* fix

* add interpolation functions

* restructure again

* fix gpu bugs

* remove unnecessary argument

* fix gpu again

* fix gpu

* implement suggestions

* fix

* add comments

* implement suggestions

* first prototype: NOT VALIDATED YET

* fix avg vel in LastiwkaModel

* fix avg in mirroring

* fix tests

* adapt example file

* fix allocations

* fix gpu

* first GPU working prototype

* improve API

* fix tests

* rm export

* fix

* fix type stability

* implement suggestions

* rm ref

* format

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* fix merge bugs

* supersede #852

* supersede #850

* fix

* rename boundary models

* fix?

* fix tests

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* fix characteristics

* revise pipe flow example

* fix include bug

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* Combine PST and TVF into a unified framework (#884)

* Combine PST and TVF into a unified framework

* Require update callback for PST

* Fix WCSPH

* Update PST only in callback

* Fix EDAC

* Update docs

* Fix alle example files

* Fix tests

* Fix periodic channel

* Fix docs

* Update news

* Fix tests

* Fix example file

* add doc strings

* implement suggestions

* remove duplicated system

* rm unnecessary variable

* apply formatter

* set pressure for WCSPH

* setter for open boundary

* rename reference values

* fix gpu tests

* fix sign

* implement suggestions

* fix unit tests

* revise doc strings

* add NEWS entry

* update NEWS

---------

Co-authored-by: LasNikas <niklas.nehe@web.de>
Co-authored-by: Sven Berger <berger.sven@gmail.com>
Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>

Author: 4 people

NEWS.md

@@ -8,6 +8,9 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 ### API Changes
 
+- API for `OpenBoundarySPHSystem` and `BoundaryZone` changed.
+  It is now possible to pass multiple `BoundaryZone`s to a single `OpenBoundarySPHSystem`.
+  Reference values are now assigned individually to each `BoundaryZone`. (#866)
 - Combined transport velocity formulation (TVF) and particle shifting technique (PST) into
   one unified framework.
   The keyword argument `transport_velocity` now changed to `shifting_technique`.

examples/fluid/pipe_flow_2d.jl

@@ -12,7 +12,7 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacing = 0.05
+particle_spacing = 0.02
 
 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled
 boundary_layers = 4
@@ -21,7 +21,7 @@ boundary_layers = 4
 # 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 = 8
+open_boundary_layers = 6
 
 # ==========================================================================================
 # ==== Experiment Setup
@@ -30,65 +30,71 @@ tspan = (0.0, 2.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 = 2.0
+const prescribed_velocity = (1.0, 0.0)
+flow_direction = [1.0, 0.0]
 
-boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
-                 domain_size[2])
+open_boundary_size = (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 = 20 * prescribed_velocity
-
-state_equation = nothing
+sound_speed = 10 * maximum(abs.(prescribed_velocity))
 
-pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
-                       pressure=pressure, n_layers=boundary_layers,
+pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density,
+                       n_layers=boundary_layers, velocity=prescribed_velocity,
                        faces=(false, false, true, true))
 
-# Shift pipe walls in negative x-direction for the inflow
-pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
+min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0)
+inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
+                        fluid_density, n_layers=boundary_layers,
+                        min_coordinates=min_coords_inlet,
+                        faces=(false, false, true, true))
+
+min_coords_outlet = (pipe.fluid_size[1], 0.0)
+outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
+                         fluid_density, n_layers=boundary_layers,
+                         min_coordinates=min_coords_outlet,
+                         faces=(false, false, true, true))
 
 NDIMS = ndims(pipe.fluid)
 
-n_buffer_particles = 5 * pipe.n_particles_per_dimension[2]^(NDIMS - 1)
+n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^(NDIMS - 1)
 
 # ==========================================================================================
 # ==== Fluid
-wcsph = false
+wcsph = true
 
 smoothing_length = 1.5 * particle_spacing
 smoothing_kernel = WendlandC2Kernel{NDIMS}()
 
 fluid_density_calculator = ContinuityDensity()
 
-kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number
+kinematic_viscosity = maximum(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,
-                                           shifting_technique=ParticleShiftingTechnique(),
-                                           buffer_size=n_buffer_particles)
-
 # Alternatively the WCSPH scheme can be used
 if wcsph
     state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
-                                       exponent=1, background_pressure=pressure)
-    alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed)
-    viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
+                                       exponent=1)
+    density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 
     fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator,
                                                state_equation, smoothing_kernel,
+                                               density_diffusion=density_diffusion,
                                                smoothing_length, viscosity=viscosity,
                                                shifting_technique=ParticleShiftingTechnique(),
                                                buffer_size=n_buffer_particles)
+else
+    # Alternatively the EDAC scheme can be used
+    state_equation = nothing
+
+    fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel,
+                                               smoothing_length, sound_speed,
+                                               viscosity=viscosity,
+                                               density_calculator=fluid_density_calculator,
+                                               shifting_technique=ParticleShiftingTechnique(),
+                                               buffer_size=n_buffer_particles)
 end
 
 # ==========================================================================================
@@ -98,63 +104,61 @@ function velocity_function2d(pos, t)
     # Use this for a time-dependent inflow velocity
     # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0)
 
-    return SVector(prescribed_velocity, 0.0)
+    return SVector(prescribed_velocity)
 end
 
-open_boundary_model = BoundaryModelLastiwka()
+open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring())
 
+reference_velocity_in = velocity_function2d
+reference_pressure_in = nothing
+reference_density_in = nothing
 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)
-
-reference_velocity_in = velocity_function2d
-reference_pressure_in = pressure
-reference_density_in = fluid_density
-open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
-                                         boundary_model=open_boundary_model,
-                                         buffer_size=n_buffer_particles,
-                                         reference_density=reference_density_in,
-                                         reference_pressure=reference_pressure_in,
-                                         reference_velocity=reference_velocity_in)
-
+                      reference_density=reference_density_in,
+                      reference_pressure=reference_pressure_in,
+                      reference_velocity=reference_velocity_in,
+                      initial_condition=inlet.fluid, boundary_type=boundary_type_in)
+
+reference_velocity_out = nothing
+reference_pressure_out = nothing
+reference_density_out = nothing
 boundary_type_out = OutFlow()
-plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]])
+plane_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[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)
-
-reference_velocity_out = velocity_function2d
-reference_pressure_out = pressure
-reference_density_out = fluid_density
-open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
-                                          boundary_model=open_boundary_model,
-                                          buffer_size=n_buffer_particles,
-                                          reference_density=reference_density_out,
-                                          reference_pressure=reference_pressure_out,
-                                          reference_velocity=reference_velocity_out)
+                       reference_density=reference_density_out,
+                       reference_pressure=reference_pressure_out,
+                       reference_velocity=reference_velocity_out,
+                       initial_condition=outlet.fluid, boundary_type=boundary_type_out)
+
+open_boundary = OpenBoundarySPHSystem(inflow, outflow; fluid_system,
+                                      boundary_model=open_boundary_model,
+                                      buffer_size=n_buffer_particles)
+
 # ==========================================================================================
 # ==== Boundary
-viscosity_boundary = ViscosityAdami(nu=1e-4)
-boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
+wall = union(pipe.boundary, inlet.boundary, outlet.boundary)
+viscosity_boundary = viscosity
+boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass,
                                              AdamiPressureExtrapolation(),
                                              state_equation=state_equation,
                                              viscosity=viscosity_boundary,
                                              smoothing_kernel, smoothing_length)
 
-boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model)
+boundary_system = BoundarySPHSystem(wall, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2)
-max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2)
+min_corner = minimum(wall.coordinates .- particle_spacing, dims=2)
+max_corner = maximum(wall.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,
+semi = Semidiscretization(fluid_system, open_boundary, boundary_system,
+                          neighborhood_search=nhs,
                           parallelization_backend=PolyesterBackend())
 
 ode = semidiscretize(semi, tspan)

examples/fluid/pipe_flow_3d.jl

@@ -22,28 +22,24 @@ open_boundary_layers = 6
 
 # ==========================================================================================
 # ==== Experiment Setup
-tspan = (0.0, 2.0)
-
-function velocity_function3d(pos, t)
-    # Use this for a time-dependent inflow velocity
-    # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0)
-
-    return SVector(prescribed_velocity, 0.0, 0.0)
-end
+tspan = (0.0, 0.5)
 
 domain_size = (1.0, 0.4, 0.4)
-
-boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
-                 domain_size[2], domain_size[3])
-
+const prescribed_velocity = (1.0, 0.0, 0.0)
 flow_direction = [1.0, 0.0, 0.0]
 
+open_boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+                      domain_size[2], domain_size[3])
+min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0)
+min_coords_outlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0)
+
 # setup simulation
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"),
-              domain_size=domain_size, boundary_size=boundary_size,
+              domain_size=domain_size, open_boundary_size=open_boundary_size,
               flow_direction=flow_direction, faces=(false, false, true, true, true, true),
-              tspan=tspan, reference_velocity=velocity_function3d,
-              open_boundary_layers=open_boundary_layers,
+              tspan=tspan, prescribed_velocity=prescribed_velocity,
+              open_boundary_layers=open_boundary_layers, min_coords_inlet=min_coords_inlet,
+              min_coords_outlet=min_coords_outlet,
               plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0],
                         [0.0, 0.0, domain_size[3]]),
               plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0],

src/TrixiParticles.jl

@@ -79,7 +79,8 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
        DensityDiffusionAntuono
 export tensile_instability_control
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
-       PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni,
+       PressureMirroring, PressureZeroing, BoundaryModelCharacteristicsLastiwka,
+       BoundaryModelMirroringTafuni,
        BernoulliPressureExtrapolation
 export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring
 export HertzContactModel, LinearContactModel

src/general/buffer.jl

@@ -36,6 +36,10 @@ function allocate_buffer(initial_condition, buffer::SystemBuffer)
     return union(initial_condition, buffer_ic)
 end
 
+# By default, there is no buffer.
+# Dispatch by system type to handle systems that provide a buffer.
+@inline buffer(system) = nothing
+
 @inline update_system_buffer!(buffer::Nothing, semi) = buffer
 
 # TODO `resize` allocates. Find a non-allocating version

src/general/semidiscretization.jl

@@ -972,17 +972,17 @@ end
 
 function check_configuration(system::OpenBoundarySPHSystem, systems,
                              neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch)
-    (; boundary_model, boundary_zone) = system
+    (; boundary_model, boundary_zones) = system
 
     # Store index of the fluid system. This is necessary for re-linking
     # in case we use Adapt.jl to create a new semidiscretization.
     fluid_system_index = findfirst(==(system.fluid_system), systems)
     system.fluid_system_index[] = fluid_system_index
 
-    if boundary_model isa BoundaryModelLastiwka &&
-       boundary_zone isa BoundaryZone{BidirectionalFlow}
-        throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " *
-                            "Please specify inflow and outflow."))
+    if boundary_model isa BoundaryModelCharacteristicsLastiwka &&
+       any(zone -> isnothing(zone.flow_direction), boundary_zones)
+        throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " *
+                            "Please specify `InFlow()` and `OutFlow()`."))
     end
 
     if first(PointNeighbors.requires_update(neighborhood_search))
@@ -1011,10 +1011,8 @@ function set_system_links(system::OpenBoundarySPHSystem, semi)
                                  system.pressure,
                                  system.boundary_candidates,
                                  system.fluid_candidates,
-                                 system.boundary_zone,
-                                 system.reference_velocity,
-                                 system.reference_pressure,
-                                 system.reference_density,
+                                 system.boundary_zone_indices,
+                                 system.boundary_zones,
                                  system.buffer,
                                  system.cache)
 end

src/general/system.jl

@@ -37,17 +37,18 @@ initialize!(system, semi) = system
 # Number of particles in the system whose positions are to be integrated (corresponds to the size of u and du)
 @inline n_moving_particles(system) = nparticles(system)
 
-@inline eachparticle(system) = Base.OneTo(nparticles(system))
+@inline eachparticle(system::System) = active_particles(system)
+@inline eachparticle(initial_condition) = Base.OneTo(nparticles(initial_condition))
 
 # Wrapper for systems with `SystemBuffer`
-@inline each_moving_particle(system) = each_moving_particle(system, system.buffer)
+@inline each_moving_particle(system) = each_moving_particle(system, buffer(system))
 @inline each_moving_particle(system, ::Nothing) = Base.OneTo(n_moving_particles(system))
 
-@inline active_coordinates(u, system) = active_coordinates(u, system, system.buffer)
+@inline active_coordinates(u, system) = active_coordinates(u, system, buffer(system))
 @inline active_coordinates(u, system, ::Nothing) = current_coordinates(u, system)
 
-@inline active_particles(system) = active_particles(system, system.buffer)
-@inline active_particles(system, ::Nothing) = eachparticle(system)
+@inline active_particles(system) = active_particles(system, buffer(system))
+@inline active_particles(system, ::Nothing) = Base.OneTo(nparticles(system))
 
 # This should not be dispatched by system type. We always expect to get a column of `A`.
 @propagate_inbounds function extract_svector(A, system, i)