Merge branch 'main' into test-1.12

Author: svchb

.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.31.2
+        uses: crate-ci/typos@v1.32.0

NEWS.md

@@ -4,6 +4,14 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.3.1
+
+### Features
+
+- With all CPU backends, a new array type is used for the integration array, which defines
+  broadcasting to be multithreaded, leading to speedups of up to 5x with large thread counts
+  when combined with thread pinning (#722).
+
 ## Version 0.3
 
 ### API Changes

Project.toml

@@ -54,7 +54,7 @@ KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
 OrdinaryDiffEq = "6.91"
 OrdinaryDiffEqCore = "1"
-PointNeighbors = "0.6.1"
+PointNeighbors = "0.6.3"
 Polyester = "0.7.10"
 ReadVTK = "0.2"
 RecipesBase = "1"

docs/src/gpu.md

@@ -13,13 +13,14 @@ Unlike the default cell list, which assumes an unbounded domain,
 this cell list requires a bounding box for the domain.
 For simulations that are bounded by a closed tank, we can simply use the boundary
 of the tank to obtain the bounding box as follows.
-```jldoctest gpu; output=false, setup=:(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing))
+```jldoctest gpu; output=false, filter = r"FullGridCellList{PointNeighbors.DynamicVectorOfVectors{.*", setup=:(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing))
 min_corner = minimum(tank.boundary.coordinates, dims=2)
 max_corner = maximum(tank.boundary.coordinates, dims=2)
 cell_list = FullGridCellList(; min_corner, max_corner)
 
 # output
-FullGridCellList{PointNeighbors.DynamicVectorOfVectors{Int32, Matrix{Int32}, Vector{Int32}, Base.RefValue{Int32}}, Nothing, SVector{2, Float64}, SVector{2, Float64}}(Vector{Int32}[], nothing, [-0.12500000000000003, -0.12500000000000003], [1.125, 1.125])
+FullGridCellList{PointNeighbors.DynamicVectorOfVectors{...}(...)
+
 ```
 
 We then need to pass this cell list to the neighborhood search and the neighborhood search

docs/src/preprocessing/preprocessing.md

@@ -1,4 +1,4 @@
-# Sampling of Geometries
+# [Sampling of Geometries](@id sampling_of_geometries)
 
 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,
@@ -250,15 +250,108 @@ Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")
 Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_jacobson.jl")]
 ```
+# [Read geometries from file](@id read_geometries_from_file)
+Geometries can be imported using the [`load_geometry`](@ref) function.
+- For 3D geometries, we support the binary (`.stl`) format.
+- For 2D geometries, the recommended format is DXF (`.dxf`), with optional support for a simple ASCII (`.asc`) format.
+
+## ASCII Format (.asc)
+An .asc file contains a list of 2D coordinates, space-delimited, one point per line,
+where the first column are the x values and the second the y values.
+For example:
+```plaintext
+# ASCII
+0.0 0.0
+1.0 0.0
+1.0 1.0
+0.0 1.0
+```
+It is the user’s responsibility to ensure the points are ordered correctly.
+This format is easy to generate and inspect manually.
+
+## DXF Format (.dxf) – recommended
+The DXF (Drawing Exchange Format) is a widely-used CAD format for 2D and 3D data.
+We recommend this format for defining 2D polygonal.
+Only DXF entities of type `LINE` or `POLYLINE` are supported.
+To create DXF files from scratch, we recommend using the open-source tool [FreeCAD](https://www.freecad.org/).
+For a less technical approach, we recommend using [Inkscape](https://inkscape.org/) to create SVG files and then export them as DXF.
+
+### Creating DXF Files with FreeCAD
+
+- Open FreeCAD and create a new document.
+- Switch to the Sketcher workbench and create a new sketch.
+- Choose the XY-plane orientation and draw your geometry.
+- Select the object to be exported.
+- Go to "File > Export..." and choose "AutoDesk DXF (*.dxf)" as the format.
+- Ensure the following Import-Export options are enabled:
+    - "Use legacy Python exporter".
+    - "Project exported objects along current view direction".
+
+### Creating DXF Files with Inkscape
+SVG vector graphics can also be used as a basis for geometry.
+Use Inkscape to open or create the SVG.
+You can simply draw a Bezier curve by pressing "B" on your keyboard.
+We reommend the mode "Create spiro paths" for a smooth curve.
+Select the relevant path and export it as DXF:
+- Go to "File > Save As...".
+- Choose "Desktop Cutting Plotter (AutoCAD DXF R12)(*.dxf)" format.
 
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 ```
 
-# Particle Packing
+# [Particle Packing](@id particle_packing)
+To obtain a body-fitted and isotropic particle distribution, an initial configuration (see [Sampling of Geometries](@ref sampling_of_geometries)) is first generated. This configuration is then packed using a [`ParticlePackingSystem`](@ref).
+The preprocessing pipeline consists of the following steps:
+
+- Load geometry: Fig. 1, [`load_geometry`](@ref).
+- Compute the signed distance field (SDF): Fig. 2, [`SignedDistanceField`](@ref).
+- Generate boundary particles: Fig. 3, [`sample_boundary`](@ref).
+- Initial sampling of the interior particles with inside-outside segmentation: Fig. 4, [`ComplexShape`](@ref).
+- Pack the initial configuration of interior and boundary particles (Fig. 5): Fig. 6, [`ParticlePackingSystem`](@ref).
+
+The input data can either be a 3D triangulated surface mesh represented in STL format or a 2D polygonal traversal of the geometry (see [`load_geometry`](@ref)).
+The second step involves generating the SDF (see [`SignedDistanceField`](@ref)), which is necessary for the final packing step as it requires a surface detection.
+The SDF is illustrated in Fig. 2, where the distances to the surface of the geometry are visualized as a color map.
+As shown, the SDF is computed only within a narrow band around the geometry’s surface, enabling  a face-based neighborhood search (NHS) to be used exclusively during this step.
+In the third step, the initial configuration of the boundary particles is generated (orange particles in Fig. 3).
+Boundary particles are created by copying the positions of SDF points located outside the geometry but within a predefined boundary thickness (see [`sample_boundary`](@ref)).
+In the fourth step, the initial configuration of the interior particles (green particles in Fig. 4) is generated using the hierarchical winding number approach (see [Hierarchical Winding](@ref hierarchical_winding)).
+After steps **1** through **4**, the initial configuration of both interior and boundary particles is obtained, as illustrated in Fig. 5.
+The interface of the geometry surface is not well resolved with the initial particle configuration.
+Thus, in the final step, a packing algorithm by Zhu et al. [Zhu2021](@cite) is applied utilizing the SDF to simultaneously optimize the positions of both interior and boundary particles,
+yielding an isotropic distribution while accurately preserving the geometry surface, as illustrated in Fig. 6.
+
+```@raw html
+<div style="display: flex; gap: 16px; flex-wrap: wrap;">
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/7fe9d1f1-1633-4377-8b97-a4d1778aee07" alt="geometry" style="max-width: 200px;">
+    <figcaption>(1) Geometry representation</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/2b79188c-3148-49f1-8337-718721851bf5" alt="sdf" style="max-width: 200px;">
+    <figcaption>(2) Signed distances to the surface</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/1501718f-d1f5-4f14-b1bc-2a2e581db476" alt="boundary" style="max-width: 200px;">
+    <figcaption>(3) Boundary particles</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/f7376b15-324a-4da1-bb59-db01c7bd6620" alt="interior" style="max-width: 200px;">
+    <figcaption>(4) Interior particles</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/4be889d6-e70a-4c5e-bef2-0071ea4d898c" alt="initial_config" style="max-width: 200px;">
+    <figcaption>(5) Initial configuration</figcaption>
+  </figure>
+  <figure style="margin: 0; text-align: center;">
+    <img src="https://github.com/user-attachments/assets/0f7aba29-3cf7-4ec1-8c95-841e72fe620d" alt="packed_config" style="max-width: 200px;">
+    <figcaption>(6) Packed configuration</figcaption>
+  </figure>
+</div>
+```
 
-TODO
 
 ```@autodocs
 Modules = [TrixiParticles]

docs/src/refs.bib

@@ -648,6 +648,17 @@ @Article{Sun2017
   publisher = {Elsevier BV},
 }
 
+@article{Sun2018,
+	author    = {Sun, P. N. and Colagrossi, A. and Marrone, S. and Antuono, M. and Zhang, A. M.},
+	title     = {Multi-resolution Delta-plus-{SPH} with tensile instability control: Towards high Reynolds number flows},
+	journal   = {Computer Physics Communications},
+	volume    = {224},
+  year      = {2018},
+	issn      = {0010-4655},
+	pages     = {63--80},
+	doi       = {10.1016/j.cpc.2017.11.016},
+}
+
 @Article{Tafuni2018,
   author    = {A. Tafuni and J.M. Dom{\'{\i}}nguez and R. Vacondio and A.J.C. Crespo},
   journal   = {Computer Methods in Applied Mechanics and Engineering},
@@ -723,4 +734,18 @@ @book{Poling2001
   year = {2001},
   publisher = {McGraw-Hill},
   address = {New York}
-}
+}
+
+@article{Zhu2021,
+  author    = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu},
+  journal   = {Journal of Hydrodynamics},
+  title     = {A {CAD}-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction},
+  year      = {2021},
+  issn      = {1878-0342},
+  month     = apr,
+  number    = {2},
+  pages     = {195--206},
+  volume    = {33},
+  doi       = {10.1007/s42241-021-0031-y},
+  publisher = {Springer Science and Business Media LLC}
+}

docs/src/systems/fluid.md

@@ -1,8 +1,8 @@
 # [Fluid Models](@id fluid_models)
 
 Currently available fluid methods are the [weakly compressible SPH method](@ref wcsph) and the
-[entropically damped artificial compressibility for SPH](@ref edac).  
-This page lists models and techniques that apply to both of these methods.  
+[entropically damped artificial compressibility for SPH](@ref edac).
+This page lists models and techniques that apply to both of these methods.
 
 ## [Viscosity](@id viscosity_sph)
 
@@ -58,21 +58,19 @@ by Balsara ([Balsara1995](@cite)) or Morris ([Morris1997](@cite)).
 The force exerted by particle `b` on particle `a` due to artificial viscosity is given by:
 
 ```math
-\mathbf{F}_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab}
+F_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab}
 ```
 
 where:
 
-- `\Pi_{ab}` is the artificial viscosity term defined as:
-
-    ```math
-    \Pi_{ab} =
-    \begin{cases}
-        -\frac{\alpha c \mu_{ab} + \beta \mu_{ab}^2}{\bar{\rho}_{ab}} & \text{if } \mathbf{v}_{ab} \cdot \mathbf{r}_{ab} < 0, \\
-        0 & \text{otherwise}
-    \end{cases}
-    ```
-
+- ``\Pi_{ab}`` is the artificial viscosity term defined as:
+  ```math
+  \Pi_{ab} =
+  \begin{cases}
+      -\frac{\alpha c \mu_{ab} + \beta \mu_{ab}^2}{\bar{\rho}_{ab}} & \text{if } \mathbf{v}_{ab} \cdot \mathbf{r}_{ab} < 0, \\
+      0 & \text{otherwise}
+  \end{cases}
+  ```
 - ``\alpha`` and ``\beta`` are viscosity parameters,
 - ``c`` is the local speed of sound,
 - ``\bar{\rho}_{ab}`` is the arithmetic mean of the densities of particles ``a`` and ``b``.
@@ -361,7 +359,7 @@ with:
 This divergence can be computed numerically in the SPH framework as
 
 ```math
-\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab} 
+\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab}
 ```
 
 #### Advantages and limitations

docs/src/systems/weakly_compressible_sph.md

@@ -54,8 +54,8 @@ pressure field. It is highly recommended to use density diffusion when using WCS
 All density diffusion terms extend the continuity equation (see [`ContinuityDensity`](@ref))
 by an additional term
 ```math
-\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla_{r_a} W(\Vert r_{ab} \Vert, h)
-    + \delta h c \sum_{b} V_b \psi_{ab} \cdot \nabla_{r_a} W(\Vert r_{ab} \Vert, h),
+\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla W_{ab}
+    + \delta h c \sum_{b} V_b \psi_{ab} \cdot \nabla W_{ab},
 ```
 where ``V_b = m_b / \rho_b`` is the volume of particle ``b`` and ``\psi_{ab}`` depends on
 the density diffusion method (see [`DensityDiffusion`](@ref) for available terms).
@@ -109,3 +109,88 @@ term is very expensive and adds about 40--50% of computational cost.
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "density_diffusion.jl")]
 ```
+
+## [Particle Shifting Technique](@id shifting)
+
+The Particle Shifting Technique (PST) is used to correct tensile instability
+in regions of low pressure, as observed in viscous flow around an object.
+Without PST, tensile instability causes non-physical separation of the fluid
+from the object, introducing void regions behind the object.
+
+At lower resolutions, PST alone can be effective to correct a viscous flow
+around a cylinder, as shown in this figure.
+![particle_shifting](https://github.com/user-attachments/assets/70892b0b-af36-4531-b328-f73e63a5e33c)
+At higher resolutions, PST alone is not effective anymore; see the figure in
+[Tensile Instability Control](@ref tic).
+We recommend using PST and Tensile Instability Control together
+in such simulations.
+
+### Mathematical formulation
+
+We use the following formulation by [Sun et al. (2018)](@cite Sun2018).
+After each time step, a correction term ``\delta r_a`` is added to the position ``r_a``
+of particle ``a``, which is given by
+```math
+\delta r_a = -4 \Delta t \, v_\text{max} h
+    \sum_b \left( 1 + R \left( \frac{W_{ab}}{W(\Delta x_a)} \right)^n \right) \nabla W_{ab}
+    \frac{m_b}{\rho_a + \rho_b},
+```
+where:
+- ``\Delta t`` is the time step,
+- ``v_\text{max}`` is the maximum velocity over all particles,
+- ``h`` is the smoothing length,
+- ``R`` and ``n`` are constants, which are set to ``0.2`` and ``4`` respectively,
+- ``W(\Delta x_a)`` is the smoothing kernel of the particle size of particle ``a``,
+  which can be interpreted as the target particle spacing that we want to achieve.
+- ``\nabla W_{ab}`` is the gradient of the smoothing kernel,
+- ``m_b`` is the mass of particle ``b``,
+- ``\rho_a, \rho_b`` is the density of particles ``a`` and ``b``, respectively.
+
+Note that we replaced ``\text{CFL} \cdot \text{Ma}`` by ``\Delta t \cdot v_\text{max} / h``,
+as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9.
+
+The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation
+of PST is commonly referred to as ``\delta^+``-SPH.
+
+The Particle Shifting Technique can be applied in form
+of the [`ParticleShiftingCallback`](@ref).
+
+## [Tensile Instability Control](@id tic)
+
+Tensile Instability Control (TIC) is a pressure acceleration formulation to correct tensile
+instability in regions of low pressure, as observed in viscous flow around an object.
+The technique was introduced by [Sun et al. (2018)](@cite Sun2018).
+The formulation is described in Section 2.1 of this paper.
+It can be used in combination with the [Particle Shifting Technique (PST)](@ref shifting)
+to effectively prevent non-physical separation of the fluid from the object.
+
+As can be seen in the following figure, TIC alone can cause instabilities
+and does not improve the simulation.
+PST alone can mostly prevent separation at lower resolutions.
+A small void region is still visible, but quickly disappears.
+The combination of PST and TIC prevents separation effectively.
+![low_res](https://github.com/user-attachments/assets/5b30d440-8ca5-4c13-94d0-d110de2eb7cc)
+
+At higher resolutions, PST alone is not effective to prevent separation, as can be seen
+in the next figure.
+Only the combination of PST and TIC is able to produce physical results.
+![high_res](https://github.com/user-attachments/assets/674aec76-33e6-4ee3-bcd7-ba8c381a2775)
+
+### Mathematical formulation
+
+The force that particle ``a`` experiences from particle ``b`` due to pressure is given by
+```math
+f_{ab} = -m_a m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla W_{ab}
+```
+for the WCSPH method with [`ContinuityDensity`](@ref).
+
+The TIC formulation changes this force to
+```math
+f_{ab} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla W_{ab}.
+```
+Note that this formulation is asymmetric and sacrifices conservation of linear and angular
+momentum.
+
+```@docs
+tensile_instability_control
+```

examples/fluid/dam_break_oil_film_2d.jl

@@ -29,7 +29,7 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 
 # TODO: broken if both systems use surface tension
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-              sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
+              sol=nothing, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan,
               viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
               gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed,
               prefix="", reference_particle_spacing=fluid_particle_spacing)

examples/fluid/periodic_channel_2d.jl

@@ -35,9 +35,11 @@ smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 
+# `pressure_acceleration=nothing` is the default and can be overwritten with `trixi_include`
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
-                                           smoothing_length, viscosity=viscosity)
+                                           smoothing_length, viscosity=viscosity,
+                                           pressure_acceleration=nothing)
 
 # ==========================================================================================
 # ==== Boundary
@@ -65,8 +67,10 @@ ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+# This can be overwritten with `trixi_include`
+extra_callback = nothing
 
-callbacks = CallbackSet(info_callback, saving_callback)
+callbacks = CallbackSet(info_callback, saving_callback, extra_callback)
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 # Limiting of the maximum stepsize is necessary to prevent crashing.