diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
index 71eea7176..fedf097ec 100644
--- a/.github/workflows/SpellCheck.yml
+++ b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
- name: Checkout Actions Repository
uses: actions/checkout@v4
- name: Check spelling
- uses: crate-ci/typos@v1.23.6
+ uses: crate-ci/typos@v1.24.3
diff --git a/README.md b/README.md
index 0f0dba7a7..1283b9c4f 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,8 @@
**TrixiParticles.jl** is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in [Julia](https://julialang.org).
A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing.
-Its features include:
+
+[](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s)
## Features
- Incompressible Navier-Stokes
diff --git a/docs/make.jl b/docs/make.jl
index 3bc2dec41..a00209908 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -104,8 +104,8 @@ makedocs(sitename="TrixiParticles.jl",
"Tutorial" => "tutorial.md",
"Examples" => "examples.md",
"Visualization" => "visualization.md",
- "preprocessing" => [
- "Sample Geometries" => joinpath("preprocessing", "preprocessing.md"),
+ "Preprocessing" => [
+ "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md"),
],
"Components" => [
"Overview" => "overview.md",
diff --git a/docs/src/install.md b/docs/src/install.md
index f88a10121..6c590f760 100644
--- a/docs/src/install.md
+++ b/docs/src/install.md
@@ -1,8 +1,8 @@
# [Installation](@id installation)
## Setting up Julia
-If you have not yet installed Julia, please [follow the instructions for your
-operating system](https://julialang.org/downloads/platform/). TrixiParticles.jl works
+If you have not yet installed Julia, please [follow the instructions on the
+official website](https://julialang.org/downloads/). TrixiParticles.jl works
with Julia v1.9 and newer. We recommend using the latest stable release of Julia.
## For users
diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md
index 8d83c9caf..67192221f 100644
--- a/docs/src/preprocessing/preprocessing.md
+++ b/docs/src/preprocessing/preprocessing.md
@@ -1,10 +1,249 @@
+# Sampling of Geometries
-TODO
+Generating the initial configuration of a simulation requires filling volumes (3D) or surfaces (2D) of complex geometries with particles.
+The algorithm to sample a complex geometry should be robust and fast,
+since for large problems (large numbers of particles) or complex geometries (many geometry faces),
+generating the initial configuration is not trivial and can be very expensive in terms of computational cost.
+We therefore use a [winding number](https://en.wikipedia.org/wiki/Winding_number) approach for an inside-outside segmentation of an object.
+The winding number ``w(\mathbf{p})`` is a signed integer-valued function of a point ``\mathbf{p}`` and is defined as
+```math
+w(\mathbf{p}) = \frac{1}{2 \pi} \sum^n_{i=1} \Theta_i.
+```
+
+Here, ``\Theta_i`` is the *signed* angle between ``\mathbf{c}_i - \mathbf{p}`` and ``\mathbf{c}_{i+1} - \mathbf{p}`` where ``\mathbf{c}_i`` and ``\mathbf{c}_{i+1}`` are two consecutive vertices on a curve.
+In 3D, we refer to the solid angle of an *oriented* triangle with respect to ``\mathbf{p}``.
+
+We provide the following methods to calculate ``w(\mathbf{p})``:
+- Hormann et al. (2001) evaluate the winding number combined with an even-odd rule, but only for 2D polygons (see [WindingNumberHormann](@ref)).
+- Naive winding: Jacobson et al. (2013) generalized the winding number so that the algorithm can be applied for both 2D and 3D geometries (see [WindingNumberJacobson](@ref)).
+- Hierarchical winding: Jacobson et al. (2013) also introduced a fast hierarchical evaluation of the winding number. For further information see the description below.
+
+## [Hierarchical Winding](@id hierarchical_winding)
+According to Jacobson et al. (2013) the winding number with respect to a polygon (2D) or triangle mesh (3D) is the sum of the winding numbers with respect to each edge (2D) or face (3D).
+We can show this with the following example in which we determine the winding number for each edge of a triangle separately and sum them up:
+
+```julia
+using TrixiParticles
+using Plots
+
+triangle = [125.0 375.0 250.0 125.0;
+ 175.0 175.0 350.0 175.0]
+
+# Delete all edges but one
+edge1 = deleteat!(TrixiParticles.Polygon(triangle), [2, 3])
+edge2 = deleteat!(TrixiParticles.Polygon(triangle), [1, 3])
+edge3 = deleteat!(TrixiParticles.Polygon(triangle), [1, 2])
+
+algorithm = WindingNumberJacobson()
+
+grid = hcat(([x, y] for x in 1:500, y in 1:500)...)
+
+_, w1 = algorithm(edge1, grid; store_winding_number=true)
+_, w2 = algorithm(edge2, grid; store_winding_number=true)
+_, w3 = algorithm(edge3, grid; store_winding_number=true)
+
+w = w1 + w2 + w3
+
+heatmap(1:500, 1:500, reshape(w1, 500, 500)', color=:coolwarm, showaxis=false,
+ tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w2, 500, 500)', color=:coolwarm, showaxis=false,
+ tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w3, 500, 500)', color=:coolwarm, showaxis=false,
+ tickfontsize=12, size=(570, 500), margin=6 * Plots.mm)
+heatmap(1:500, 1:500, reshape(w, 500, 500)', color=:coolwarm, showaxis=false,
+ tickfontsize=12, size=(570, 500), margin=6 * Plots.mm, clims=(-1, 1))
+
+```
+
+```@raw html
+
+
+
+```
+
+This summation property has some interesting consequences that we can utilize for an efficient computation of the winding number.
+Let ``\mathcal{S}`` be an open surface and ``\bar{\mathcal{S}}`` an arbitrary closing surface, such that
+
+```math
+\partial \bar{\mathcal{S}} = \partial \mathcal{S}
+```
+
+and ``\mathcal{B} = \bar{\mathcal{S}} \cup \mathcal{S}`` is some closed oriented surface.
+For any query point ``\mathbf{p}`` outside of ``\mathcal{B}``, we know that
+
+```math
+w_{\mathcal{S}}(\mathbf{p}) + w_{\bar{\mathcal{S}}}(\mathbf{p}) = w_{\mathcal{B}}(\mathbf{p}) = 0.
+```
+
+This means
+
+```math
+w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p}),
+```
+
+regardless of how ``\bar{\mathcal{S}}`` is constructed (as long as ``\mathbf{p}`` is outside of ``\mathcal{B}``).
+
+We can use this property in the discrete case to efficiently compute the winding number of a query point
+by partitioning the polygon or mesh in a "small" part (as in consisting of a small number of edges/faces) and a "large" part.
+For the small part we just compute the winding number, and for the large part we construct a small closing and compute its winding number.
+The partitioning is based on a hierarchical construction of bounding boxes.
+
+### Bounding volume hierarchy
+
+To efficiently find a "small part" and a "large part" as mentioned above, we construct a hierarchy of bounding boxes by starting with the whole domain and recursively splitting it in two equally sized boxes.
+The resulting hierarchy is a binary tree.
+
+The algorithm by Jacobsen et al. (Algorithm 2, p. 5) traverses this binary tree recursively until we find the leaf in which the query point is located.
+The recursion stops with the following criteria:
+
+- if the bounding box ``T`` is a leaf then ``T.\mathcal{S} = \mathcal{S} \cap T``, the part of ``\mathcal{S}``
+ that lies inside ``T``, is the "small part" mentioned above, so evaluate the winding number naively as ``w(\mathbf{p}, T.\mathcal{S})``.
+- else if ``\mathbf{p}`` is outside ``T`` then ``T.\mathcal{S}`` is the "large part", so evaluate the winding number naively
+ as ``-w(\mathbf{p}, T.\bar{\mathcal{S}})``, where ``T.\bar{\mathcal{S}}`` is the closing surface of ``T.\mathcal{S}``.
+
+#### Continuous example
+
+Now consider the following continuous (not discretized to a polygon) 2D example.
+We compute the winding number of the point ``\mathbf{p}`` with respect to ``\mathcal{S}`` using the depicted hierarchy of bounding boxes.
+
+```@raw html
+
+
+
+```
+
+(1):
+- Recurse left: ``w_{\text{left}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p}, T.\text{left})``
+- Recurse right: ``w_{\text{right}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p},T.\text{right})``
+
+(2):
+- Query point ``\mathbf{p}`` is outside bounding box ``T``, so don't recurse deeper.
+- Compute ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})`` with the closure ``T.\bar{\mathcal{S}}``, which is generally much smaller (fewer edges in the discrete version) than ``T.\mathcal{S}``:
+
+```math
+w_{\text{left}} = -\text{\texttt{naive\_winding}} (\mathbf{p}, T.\bar{\mathcal{S}})
+```
+
+(3):
+- Bounding box ``T`` is a leaf. Use open surface ``T.\mathcal{S}``:
+
+```math
+w_{\text{right}} = \text{\texttt{naive\_winding}} (\mathbf{p}, T.\mathcal{S})
+```
+
+The reconstructed surface will then look as in the following image.
+
+```@raw html
+
+
+
+```
+
+We finally sum up the winding numbers
+
+```math
+w = w_{\text{left}} + w_{\text{right} } = -w_{T_{\text{left}}.\bar{\mathcal{S}}} + w_{T_{\text{right}}.\mathcal{S}}
+```
+
+#### Discrete example
+
+We will now go through the discrete version of the example above.
+
+```@raw html
+
+
+
+```
+
+To construct the hierarchy for the discrete piecewise-linear example in (1), we have to do the following.
+
+(2):
+Each edge is distributed to the child whose box contains the edge's barycenter (red dots in (2)).
+Splitting stops when the number of a box's edges slips below a
+threshold (usually ``\approx 100`` faces in 3D, here: 6 edges).
+
+(3):
+For the closure, Jacobson et al. (2013) define *exterior vertices* (*exterior edges* in 3D)
+as boundary vertices of such a segmentation (red dots in (3)).
+To find them, we traverse around each edge (face in 3D) in order, and
+increment or decrement for each vertex (edge) a specific counter.
+
+```julia
+v1 = edge_vertices_ids[edge][1]
+v2 = edge_vertices_ids[edge][2]
+
+vertex_count[v1] += 1
+vertex_count[v2] -= 1
+```
+
+In 2D, a vertex is declared as exterior if `vertex_count(vertex) != 0`, so there is not the same amount of edges in this box going into versus out of the vertex.
+To construct the closing surface, the exterior vertices are then connected to one arbitrary
+exterior vertex using appropriately oriented line segments:
+
+```julia
+edge = vertex_count[v] > 0 ? (closing_vertex, v) : (v, closing_vertex)
+```
+
+The resulting closed surface ``T.S \cup T.\bar{S}`` then has the same number of edges going into and out of each vertex.
+
+#### Incorrect evaluation
+
+If we follow the algorithm, we know that recursion stops if
+
+- the bounding box ``T`` is a leaf or
+- the query point ``\mathbf{p}`` is outside the box.
+
+```@raw html
+
+
+
+```
+
+(1): The query point ``\mathbf{p}`` is outside the box, so we calculate the winding number with the (red) closure of the box.
+
+(2): The query point ``\mathbf{p}`` is inside the box, so we use the (blue) edges distributed to the box.
+
+(3): In this case, it leads to an incorrect evaluation of the winding number.
+The query point is clearly inside the box, but not inside the reconstructed surface.
+This is because the property ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})``
+only holds when ``\mathbf{p}`` is outside of ``\mathcal{B}``, which is not the case here.
+
+#### Correct evaluation
+Jacobson et al. (2013) don't mention this problem or provide a solution to it.
+We contacted the authors and found that they know about this problem and solve it
+by resizing the bounding box to fully include the closing surface
+of the neighboring box, since it doesn't matter if the boxes overlap.
+
+```@raw html
+
+
+
+```
+
+To avoid resizing, we take a different approach and calculate the closure of the bounding box differently:
+- Exclude intersecting edges in the calculation of the exterior vertices.
+- This way, all exterior vertices are inside the bounding box, and so will be the closing surface.
+- The intersecting edges are later added with flipped orientation,
+ so that the closing is actually a closing of the exterior plus intersecting edges.
+
+```@raw html
+
+
+
+```
+
+The evaluation then looks as follows.
+
+```@raw html
+
+
+
+```
```@autodocs
Modules = [TrixiParticles]
-Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_horman.jl")]
+Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")]
```
```@autodocs
@@ -17,11 +256,10 @@ Modules = [TrixiParticles]
Pages = [joinpath("preprocessing", "geometries", "io.jl")]
```
-
### [References](@id references_complex_shape)
- Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung "Robust inside-outside segmentation using generalized winding numbers".
In: ACM Transactions on Graphics, 32.4 (2013), pages 1--12.
[doi: 10.1145/2461912.2461916](https://igl.ethz.ch/projects/winding-number/robust-inside-outside-segmentation-using-generalized-winding-numbers-siggraph-2013-jacobson-et-al.pdf)
-- Kai Horman, Alexander Agathos "The point in polygon problem for arbitrary polygons".
+- Kai Hormann, Alexander Agathos "The point in polygon problem for arbitrary polygons".
In: Computational Geometry, 20.3 (2001), pages 131--144.
[doi: 10.1016/s0925-7721(01)00012-8](https://doi.org/10.1016/S0925-7721(01)00012-8)
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/TrixiParticles.jl b/src/TrixiParticles.jl
index 637a8bc08..d6e0c2996 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -73,7 +73,7 @@ export BoundaryMovement
export examples_dir, validation_dir, trixi_include
export trixi2vtk
export RectangularTank, RectangularShape, SphereShape, ComplexShape
-export WindingNumberHorman, WindingNumberJacobson
+export WindingNumberHormann, WindingNumberJacobson
export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry
export SourceTermDamping
export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
diff --git a/src/general/buffer.jl b/src/general/buffer.jl
index b12cf05f3..aae5ccb5c 100644
--- a/src/general/buffer.jl
+++ b/src/general/buffer.jl
@@ -1,10 +1,11 @@
struct SystemBuffer{V}
- active_particle :: BitVector
+ active_particle :: Vector{Bool}
eachparticle :: V # Vector{Int}
buffer_size :: Int
function SystemBuffer(active_size, buffer_size::Integer)
- active_particle = vcat(trues(active_size), falses(buffer_size))
+ # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
+ active_particle = vcat(fill(true, active_size), fill(false, buffer_size))
eachparticle = collect(1:active_size)
return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size)
diff --git a/src/preprocessing/geometries/triangle_mesh.jl b/src/preprocessing/geometries/triangle_mesh.jl
index e8cc993af..60d850e81 100644
--- a/src/preprocessing/geometries/triangle_mesh.jl
+++ b/src/preprocessing/geometries/triangle_mesh.jl
@@ -194,7 +194,8 @@ function unique_sorted(vertices)
# Sort by the first entry of the vectors
compare_first_element = (x, y) -> x[1] < y[1]
vertices_sorted = sort!(vertices, lt=compare_first_element)
- keep = trues(length(vertices_sorted))
+ # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
+ keep = fill(true, length(vertices_sorted))
PointNeighbors.@threaded vertices_sorted for i in eachindex(vertices_sorted)
# We only sorted by the first entry, so we have to check all previous vertices
diff --git a/src/preprocessing/point_in_poly/point_in_poly.jl b/src/preprocessing/point_in_poly/point_in_poly.jl
index 8f799757e..70ebaaf90 100644
--- a/src/preprocessing/point_in_poly/point_in_poly.jl
+++ b/src/preprocessing/point_in_poly/point_in_poly.jl
@@ -1,3 +1,3 @@
include("hierarchical_winding.jl")
-include("winding_number_horman.jl")
+include("winding_number_hormann.jl")
include("winding_number_jacobson.jl")
diff --git a/src/preprocessing/point_in_poly/winding_number_horman.jl b/src/preprocessing/point_in_poly/winding_number_hormann.jl
similarity index 79%
rename from src/preprocessing/point_in_poly/winding_number_horman.jl
rename to src/preprocessing/point_in_poly/winding_number_hormann.jl
index b1b6f9cba..02e835cd7 100644
--- a/src/preprocessing/point_in_poly/winding_number_horman.jl
+++ b/src/preprocessing/point_in_poly/winding_number_hormann.jl
@@ -1,20 +1,18 @@
"""
- WindingNumberHorman()
+ WindingNumberHormann()
-Algorithm for inside-outside segmentation of a complex geometry proposed by [Horman et al. (2001)](@ref references_complex_shape).
+Algorithm for inside-outside segmentation of a complex geometry proposed by [Hormann et al. (2001)](@ref references_complex_shape).
It is only supported for 2D geometries.
-[`WindingNumberHorman`](@ref) might handle edge cases a bit better, since the winding number is an integer value.
-Also, it is faster than [`WindingNumberJacobson`](@ref) for 2D geometries with more than about 100 edges.
-
+[`WindingNumberHormann`](@ref) might handle edge cases a bit better, since the winding number is an integer value.
!!! warning "Experimental Implementation"
This is an experimental feature and may change in any future releases.
"""
-struct WindingNumberHorman end
+struct WindingNumberHormann end
-# Algorithm 2 from Horman et al. (2001) "The point in polygon problem for arbitrary polygons"
+# Algorithm 2 from Hormann et al. (2001) "The point in polygon problem for arbitrary polygons"
# https://doi.org/10.1016/S0925-7721(01)00012-8
-function (point_in_poly::WindingNumberHorman)(geometry, points; store_winding_number=false)
+function (point_in_poly::WindingNumberHormann)(geometry, points; store_winding_number=false)
(; edge_vertices) = geometry
# We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
diff --git a/src/preprocessing/point_in_poly/winding_number_jacobson.jl b/src/preprocessing/point_in_poly/winding_number_jacobson.jl
index 888df222c..fabf2f541 100644
--- a/src/preprocessing/point_in_poly/winding_number_jacobson.jl
+++ b/src/preprocessing/point_in_poly/winding_number_jacobson.jl
@@ -58,7 +58,7 @@ Algorithm for inside-outside segmentation of a complex geometry proposed by Jaco
- `geometry`: Complex geometry returned by [`load_geometry`](@ref) and is only required when using
`hierarchical_winding=true`.
- `hierarchical_winding`: If set to `true`, an optimized hierarchical approach will be used,
- which gives a significant speedup.
+ which gives a significant speedup. For further information see [Hierarchical Winding](@ref hierarchical_winding).
- `winding_number_factor`: For leaky geometries, a factor of `0.4` will give a better inside-outside segmentation.
!!! warning "Experimental Implementation"
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 586d565b1..d8ec0a051 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -326,7 +326,7 @@ end
function remove_outside_particles(initial_condition, spanning_set, zone_origin)
(; coordinates, density, particle_spacing) = initial_condition
- in_zone = trues(nparticles(initial_condition))
+ in_zone = fill(true, nparticles(initial_condition))
for particle in eachparticle(initial_condition)
current_position = current_coords(coordinates, initial_condition, particle)
diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl
index 39613f4b1..38fdd477a 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.fluid_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.fluid_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..56e7746fc 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
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,
diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl
index fb5dcb1e7..4b1592ce6 100644
--- a/src/setups/complex_shape.jl
+++ b/src/setups/complex_shape.jl
@@ -10,7 +10,7 @@
Sample a complex geometry with particles. Returns an [`InitialCondition`](@ref).
Note that an initial particle grid is generated inside the bounding box of the geometry.
A `point_in_geometry_algorithm` checks if particles are inside the geometry or not.
-For more information about the method see [`WindingNumberJacobson`](@ref) or [`WindingNumberHorman`](@ref).
+For more information about the method see [`WindingNumberJacobson`](@ref) or [`WindingNumberHormann`](@ref).
# Arguments
- `geometry`: Geometry returned by [`load_geometry`](@ref).
@@ -30,7 +30,7 @@ For more information about the method see [`WindingNumberJacobson`](@ref) or [`W
will be overwritten when using an initial pressure function in the system.
- `point_in_geometry_algorithm`: Algorithm for sampling the complex geometry with particles.
It basically checks whether a particle is inside an object or not.
- For more information see [`WindingNumberJacobson`](@ref) or [`WindingNumberHorman`](@ref)
+ For more information see [`WindingNumberJacobson`](@ref) or [`WindingNumberHormann`](@ref)
- `grid_offset`: Offset of the initial particle grid of the bounding box of the `geometry`.
- `max_nparticles`: Maximum number of particles in the initial particle grid.
This is only used to avoid accidentally choosing a `particle_spacing`
@@ -48,8 +48,8 @@ function ComplexShape(geometry::Union{TriangleMesh, Polygon}; particle_spacing,
winding_number_factor=sqrt(eps())),
store_winding_number=false, grid_offset::Real=0.0,
max_nparticles=10^7, pad_initial_particle_grid=2particle_spacing)
- if ndims(geometry) == 3 && point_in_geometry_algorithm isa WindingNumberHorman
- throw(ArgumentError("`WindingNumberHorman` only supports 2D geometries"))
+ if ndims(geometry) == 3 && point_in_geometry_algorithm isa WindingNumberHormann
+ throw(ArgumentError("`WindingNumberHormann` only supports 2D geometries"))
end
if grid_offset < 0.0
diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl
index dabd463dc..63d7628f4 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_)
@@ -325,7 +330,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)
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(),
diff --git a/test/preprocessing/data/coordinates_horman_circle.csv b/test/preprocessing/data/coordinates_hormann_circle.csv
similarity index 100%
rename from test/preprocessing/data/coordinates_horman_circle.csv
rename to test/preprocessing/data/coordinates_hormann_circle.csv
diff --git a/test/preprocessing/data/coordinates_horman_hexagon.csv b/test/preprocessing/data/coordinates_hormann_hexagon.csv
similarity index 100%
rename from test/preprocessing/data/coordinates_horman_hexagon.csv
rename to test/preprocessing/data/coordinates_hormann_hexagon.csv
diff --git a/test/preprocessing/data/coordinates_horman_inverted_open_curve.csv b/test/preprocessing/data/coordinates_hormann_inverted_open_curve.csv
similarity index 100%
rename from test/preprocessing/data/coordinates_horman_inverted_open_curve.csv
rename to test/preprocessing/data/coordinates_hormann_inverted_open_curve.csv
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,
diff --git a/test/setups/complex_shape.jl b/test/setups/complex_shape.jl
index 977707bc3..de53aea6f 100644
--- a/test/setups/complex_shape.jl
+++ b/test/setups/complex_shape.jl
@@ -4,7 +4,7 @@
@testset verbose=true "2D" begin
@testset verbose=true "Shifted Rectangle" begin
- algorithms = [WindingNumberHorman(), WindingNumberJacobson()]
+ algorithms = [WindingNumberHormann(), WindingNumberJacobson()]
shifts = [-0.5, 0.0, 0.5]
particle_spacings = [0.03, 0.05]
@@ -39,8 +39,8 @@
@testset verbose=true "Real World Data" begin
files = ["hexagon", "circle", "inverted_open_curve"]
- algorithms = [WindingNumberHorman(), WindingNumberJacobson()]
- algorithm_names = ["horman", "jacobson"]
+ algorithms = [WindingNumberHormann(), WindingNumberJacobson()]
+ algorithm_names = ["hormann", "jacobson"]
@testset verbose=true "Algorithm: $(TrixiParticles.type2string(algorithms[i]))" for i in 1:2
@testset verbose=true "Test File `$(files[j])`" for j in eachindex(files)
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,