Merge remote-tracking branch 'refs/remotes/origin/morris_surface_tens…

…ion'

Conflicts:
	examples/fluid/waterfall.jl
	src/schemes/boundary/open_boundary/method_of_characteristics.jl
	src/schemes/boundary/open_boundary/system.jl
	src/visualization/write2vtk.jl

.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.25.0

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:
+
+[![YouTube](https://github.com/user-attachments/assets/dc2be627-a799-4bfd-9226-2077f737c4b0)](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s)
 
 ## Features
 - Incompressible Navier-Stokes

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

docs/src/preprocessing/preprocessing.md

@@ -15,7 +15,7 @@ Here, ``\Theta_i`` is the *signed* angle between ``\mathbf{c}_i - \mathbf{p}`` a
 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})``:
-- Horman et al. (2001) evaluate the winding number combined with an even-odd rule, but only for 2D polygons (see [WindingNumberHorman](@ref)).
+- 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.
 
@@ -243,7 +243,7 @@ The evaluation then looks as follows.
 
 ```@autodocs
 Modules = [TrixiParticles]
-Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_horman.jl")]
+Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")]
 ```
 
 ```@autodocs
@@ -260,6 +260,6 @@ Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 - 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)

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

examples/fluid/dam_break_oil_film_2d.jl

@@ -25,14 +25,22 @@ nu_ratio = nu_water / nu_oil
 nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
 nu_sim_water = nu_ratio * nu_sim_oil
 
+oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
+# Physical values
+nu_water = 8.9E-7
+nu_oil = 6E-5
+nu_ratio = nu_water / nu_oil
+
+nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
+nu_sim_water = nu_ratio * nu_sim_oil
+
 oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 
 # TODO: broken if both are set to surface tension
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
               viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
-              gravity=gravity, tspan=tspan,
-              density_diffusion=nothing,
+              gravity=gravity, tspan=tspan, density_diffusion=nothing,
               sound_speed=sound_speed, prefix="")
 
 # ==========================================================================================

examples/fluid/waterfall.jl

@@ -96,7 +96,7 @@ boundary_system2 = BoundarySPHSystem(end_tank.boundary, boundary_model2)
 
 outflow = OutFlow(; plane=([0.4, -0.2, -0.1], [0.4, -0.2, 0.3], [0.8, -0.2, 0.3]),
                   flow_direction=[0.0, -1.0, 0.0], open_boundary_layers=1,
-                  density=2*eps(), particle_spacing=fluid_particle_spacing)
+                  density=2 * eps(), particle_spacing=fluid_particle_spacing)
 
 open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
                                           boundary_model=BasicOutlet())

src/TrixiParticles.jl

@@ -74,7 +74,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,

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)

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

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

src/preprocessing/point_in_poly/winding_number_horman.jl → src/preprocessing/point_in_poly/winding_number_hormann.jl

@@ -1,18 +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.
+[`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

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -239,7 +239,7 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD}
         elseif !isnothing(extrude_geometry) && open_boundary_layers > 0
             initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
                                                                 particle_spacing, density,
-                                                                direction=-flow_direction_,
+                                                                direction=flow_direction_,
                                                                 n_extrude=open_boundary_layers)
         else
             initial_condition = InitialCondition(; coordinates=zeros(ELTYPE, NDIMS, 0), density=zero(ELTYPE), particle_spacing=zero(ELTYPE))
@@ -346,7 +346,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)

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)
@@ -133,6 +135,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,

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.

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,10 +195,17 @@ 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::BoundaryDEMSystem) = system.coordinates
+
+@inline initial_coordinates(system::BoundarySPHSystem, ::Nothing) = system.coordinates
+
+# We need static initial coordinates as reference when system is moving
+@inline function initial_coordinates(system::BoundarySPHSystem, movement)
+    return system.cache.initial_coordinates
 end
 
 function (movement::BoundaryMovement)(system, t)

src/schemes/fluid/surface_normal_sph.jl

@@ -57,8 +57,7 @@ function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, u_syste
                                             neighbor_system, neighbor)
         grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length)
         for i in 1:ndims(system)
-            cache.surface_normal[i, particle] += m_b / density_neighbor *
-                                                 grad_kernel[i]
+            cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i]
         end
 
         cache.neighbor_count[particle] += 1
@@ -79,17 +78,17 @@ function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_sy
     neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system)
     nhs = get_neighborhood_search(system, neighbor_system, semi)
 
-    # if smoothing_length != system.smoothing_length ||
-    #    smoothing_kernel !== system.smoothing_kernel
-    #     # TODO: this is really slow but there is no way to easily implement multiple search radia
-    #     search_radius = compact_support(smoothing_kernel, smoothing_length)
-    #     nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius,
-    #                                                   nparticles(system))
-    #     nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius,
-    #                                                   nparticles(neighbor_system))
-    #     PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords)
-    #     PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords)
-    # end
+    if smoothing_length != system.smoothing_length ||
+       smoothing_kernel !== system.smoothing_kernel
+        # TODO: this is really slow but there is no way to easily implement multiple search radia
+        search_radius = compact_support(smoothing_kernel, smoothing_length)
+        nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius,
+                                                      nparticles(system))
+        nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius,
+                                                          nparticles(neighbor_system))
+        PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords)
+        PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords)
+    end
 
     # First we need to calculate the smoothed colorfield values
     # TODO: move colorfield to extra step
@@ -118,8 +117,7 @@ function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_sy
             grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance,
                                       smoothing_length)
             for i in 1:ndims(system)
-                cache.surface_normal[i, particle] += m_b / density_neighbor *
-                                                     grad_kernel[i]
+                cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i]
             end
             cache.neighbor_count[particle] += 1
         end

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
 
@@ -244,18 +238,16 @@ 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 * system_sound_speed(particle_system) *
                          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