Adaptive Particle Refinement: Split particles

.github/workflows/Documenter.yml

@@ -25,7 +25,10 @@ concurrency:
 jobs:
   build-docs:
     name: Build and Deploy Documentation
-    runs-on: ubuntu-latest
+    runs-on: ${{ matrix.os }}
+    strategy:
+            matrix:
+                os: [ubuntu-latest, windows-latest]
     steps:
       - name: Check out project
         uses: actions/checkout@v4
@@ -40,7 +43,11 @@ jobs:
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
+        if: matrix.os == 'ubuntu-latest'
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
         run: julia --project=docs --color=yes docs/make.jl
+      - name: Just Build
+        if: matrix.os != 'ubuntu-latest'
+        run: julia --project=docs --color=yes docs/make.jl

README.md

@@ -10,8 +10,14 @@
 [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10797541.svg)](https://zenodo.org/doi/10.5281/zenodo.10797541)
 
-**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.
+**TrixiParticles.jl** is a high-performance numerical simulation framework for particle-based methods, focused on the simulation of complex multiphysics problems, and written in [Julia](https://julialang.org).
+
+TrixiParticles.jl focuses on the following use cases:
+- Accurate and efficient physics-based modelling of complex multiphysics problems.
+- Development of new particle-based methods and models.
+- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work.
+
+It offers intuitive configuration, robust pre- and post-processing, and vendor-agnostic GPU-support based on the Julia package [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). 
 
 [![YouTube](https://github.com/user-attachments/assets/dc2be627-a799-4bfd-9226-2077f737c4b0)](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s)
 

docs/make.jl

@@ -19,9 +19,12 @@ function copy_file(filename, replaces...;
     content = read(source_path, String)
     content = replace(content, replaces...)
 
+    # Use `replace` to make sure the path uses forward slashes for URLs
+    filename_url = replace(filename, "\\" => "/")
+
     header = """
     ```@meta
-    EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename"
+    EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename_url"
     ```
     """
     content = header * content
@@ -117,7 +120,8 @@ makedocs(sitename="TrixiParticles.jl",
              "Preprocessing" => [
                  "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md")
              ],
-             "Components" => [
+             "GPU Support" => "gpu.md",
+             "API Reference" => [
                  "Overview" => "overview.md",
                  "General" => [
                      "Semidiscretization" => joinpath("general", "semidiscretization.md"),

docs/src/gpu.md

@@ -0,0 +1,61 @@
+# GPU Support
+
+GPU support is still an experimental feature that is actively being worked on.
+As of now, the [`WeaklyCompressibleSPHSystem`](@ref) and the [`BoundarySPHSystem`](@ref)
+are supported on GPUs.
+We have tested this on GPUs by Nvidia and AMD.
+
+To run a simulation on a GPU, we need to use the [`FullGridCellList`](@ref)
+as cell list for the [`GridNeighborhoodSearch`](@ref).
+This cell list requires a bounding box for the domain, unlike the default cell list, which
+uses an unbounded domain.
+For simulations that are bounded by a closed tank, we can 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))
+search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length)
+min_corner = minimum(tank.boundary.coordinates, dims=2) .- search_radius
+max_corner = maximum(tank.boundary.coordinates, dims=2) .+ search_radius
+cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, max_corner)
+
+# output
+PointNeighbors.FullGridCellList{PointNeighbors.DynamicVectorOfVectors{Int32, Matrix{Int32}, Vector{Int32}, Base.RefValue{Int32}}, Nothing, SVector{2, Float64}, SVector{2, Float64}}(Vector{Int32}[], nothing, [-0.24500000000000002, -0.24500000000000002], [1.245, 1.245])
+```
+
+We then need to pass this cell list to the neighborhood search and the neighborhood search
+to the [`Semidiscretization`](@ref).
+```jldoctest gpu; output=false
+semi = Semidiscretization(fluid_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch{2}(; cell_list))
+
+# output
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ Semidiscretization                                                                               │
+│ ══════════════════                                                                               │
+│ #spatial dimensions: ………………………… 2                                                                │
+│ #systems: ……………………………………………………… 2                                                                │
+│ neighborhood search: ………………………… GridNeighborhoodSearch                                           │
+│ total #particles: ………………………………… 636                                                              │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+```
+
+At this point, we should run the simulation and make sure that it still works and that
+the bounding box is large enough.
+For some simulations where particles move outside the initial tank coordinates,
+for example when the tank is not closed or when the tank is moving, an appropriate
+bounding box has to be specified.
+
+Then, we only need to specify the data type that is used for the simulation.
+On an Nvidia GPU, we specify:
+```julia
+using CUDA
+ode = semidiscretize(semi, tspan, data_type=CuArray)
+```
+On an AMD GPU, we use:
+```julia
+using AMDGPU
+ode = semidiscretize(semi, tspan, data_type=ROCArray)
+```
+Then, we can run the simulation as usual.
+All data is transferred to the GPU during initialization and all loops over particles
+and their neighbors will be executed on the GPU as kernels generated by KernelAbstractions.jl.
+Data is only copied to the CPU for saving VTK files via the [`SolutionSavingCallback`](@ref).

docs/src/index.md

@@ -1,6 +1,14 @@
 # TrixiParticles.jl
 
-TrixiParticles.jl is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in Julia. 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:
+**TrixiParticles.jl** is a high-performance particle simulation framework designed to overcome challenges of particle-based numerical methods in multiphysics applications. Existing frameworks often lack user-friendliness, involve complex configuration, and are not easily extensible for development of new methods. In the future we also want to provide seamless scalability from CPU to Exascale-level computing with GPU support. **TrixiParticles.jl** addresses these limitations with an intuitive interface, straightforward configuration, and an extensible design, facilitating efficient simulation setup and execution.
+
+TrixiParticles.jl focuses on the following use cases:
+
+- Development of new particle-based methods and models. By providing an extensible architecture to incorporate additional particle methods easily and not focusing on a single model or numerical method.
+- Accurate, reliable and efficient physics-based modelling of complex multiphysics problems by providing a flexible configuration system, tools, high performance and a wide range of validation and test cases.
+- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work through extensive documentation, community engagement and readable configuration files.
+
+Its features include:
 
 ## Features
 - Incompressible Navier-Stokes

docs/src/overview.md

@@ -1,10 +1,16 @@
 # Overview
+The actual API reference is not listed on a single page, like in most Julia packages,
+but instead is split into multiple sections that follow a similar structure
+as the code files themselves.
+In these sections, API docs are combined with explanations of the theoretical background
+of these methods.
+
 The following page gives a rough overview of important parts of the code.
 
 ## Program flow
 
 To initiate a simulation, the goal is to solve an ordinary differential equation, for example,
-by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then 
+by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then
 utilized to integrate ``\mathrm{d}u/\mathrm{d}t`` and ``\mathrm{d}v/\mathrm{d}t``, where ``u``
 represents the particles' positions and ``v`` their properties such as velocity and density.
 During a single time step or an intermediate step of the time integration scheme, the functions

docs/src/reference-pointneighbors.md

@@ -1,4 +1,4 @@
-# PointNeighbors.jl API
+# [PointNeighbors.jl API](@id pointneighbors)
 
 ```@meta
 CurrentModule = PointNeighbors

docs/src/systems/boundary.md

@@ -55,24 +55,28 @@ of the boundary particle ``b``.
 
 ### Hydrodynamic density of dummy particles
 
-We provide five options to compute the boundary density and pressure, determined by the `density_calculator`:
+We provide six options to compute the boundary density and pressure, determined by the `density_calculator`:
 1. (Recommended) With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the
    fluid according to [Adami et al., 2012](@cite Adami2012), and the density is obtained by applying the inverse of the state equation.
    This option usually yields the best results of the options listed here.
-2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
+2. (Only relevant for FSI) With [`BernoulliPressureExtrapolation`](@ref), the pressure is extrapolated from the 
+   pressure similar to the [`AdamiPressureExtrapolation`](@ref), but a relative velocity-dependent pressure part
+   is calculated between moving solids and fluids, which increases the boundary pressure in areas prone to 
+   penetrations.
+3. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
    and the pressure is computed from the density with the state equation.
-3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
+4. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
    and the pressure is computed from the density with the state equation.
    Note that this causes a gap between fluid and boundary where the boundary is initialized
    without any contact to the fluid. This is due to overestimation of the boundary density
    as soon as the fluid comes in contact with boundary particles that initially did not have
    contact to the fluid.
    Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat.
    See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13).
-4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure
+5. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure
    is computed from the density with the state equation.
    This option is not recommended. The other options yield significantly better results.
-5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure
+6. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure
    is not used. Instead, the fluid pressure is mirrored as boundary pressure in the
    momentum equation.
    This option is not recommended due to stability issues. See [`PressureMirroring`](@ref)
@@ -93,7 +97,20 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den
     AdamiPressureExtrapolation
 ```
 
-#### 4. [`PressureZeroing`](@ref)
+#### 2. [`BernoulliPressureExtrapolation`](@ref)
+Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation:
+```math
+p_b = \frac{\sum_f (p_f + \frac{1}{2} \, \rho_{\text{neighbor}} \left( \frac{ (\mathbf{v}_f - \mathbf{v}_{\text{body}}) \cdot (\mathbf{x}_f - \mathbf{x}_{\text{neighbor}}) }{ \left\| \mathbf{x}_f - \mathbf{x}_{\text{neighbor}} \right\| } \right)^2 \times \text{factor} +\rho_f (\bm{g} - \bm{a}_b) \cdot \bm{r}_{bf}) W(\Vert r_{bf} \Vert, h)}{\sum_f W(\Vert r_{bf} \Vert, h)} 
+```
+where ``\mathbf{v}_f`` is the velocity of the fluid and ``\mathbf{v}_{\text{body}}`` is the velocity of the body.
+This adjustment provides a higher boundary pressure for solid bodies moving with a relative velocity to the fluid to prevent penetration.
+This modification is original and not derived from any literature source.
+
+```@docs
+    BernoulliPressureExtrapolation
+```
+
+#### 5. [`PressureZeroing`](@ref)
 
 This is the simplest way to implement dummy boundary particles.
 The density of each particle is set to the reference density and the pressure to the
@@ -102,7 +119,7 @@ reference pressure (the corresponding pressure to the reference density by the s
     PressureZeroing
 ```
 
-#### 5. [`PressureMirroring`](@ref)
+#### 6. [`PressureMirroring`](@ref)
 
 Instead of calculating density and pressure for each boundary particle, we modify the
 momentum equation,

examples/fluid/pipe_flow_2d.jl

@@ -48,7 +48,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi
 # Shift pipe walls in negative x-direction for the inflow
 pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
 
-n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]
+n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1)
 
 # ==========================================================================================
 # ==== Fluid
@@ -77,24 +77,27 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi
 
 # ==========================================================================================
 # ==== Open Boundary
-function velocity_function(pos, t)
+
+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)
 end
 
-inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
+plane_in = ([0.0, 0.0], [0.0, domain_size[2]])
+inflow = InFlow(; plane=plane_in, flow_direction,
                 open_boundary_layers, density=fluid_density, particle_spacing)
 
 open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
                                          boundary_model=BoundaryModelLastiwka(),
                                          buffer_size=n_buffer_particles,
                                          reference_density=fluid_density,
                                          reference_pressure=pressure,
-                                         reference_velocity=velocity_function)
+                                         reference_velocity=velocity_function2d)
 
-outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
+plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]])
+outflow = OutFlow(; plane=plane_out,
                   flow_direction, open_boundary_layers, density=fluid_density,
                   particle_spacing)
 
@@ -103,7 +106,7 @@ open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
                                           buffer_size=n_buffer_particles,
                                           reference_density=fluid_density,
                                           reference_pressure=pressure,
-                                          reference_velocity=velocity_function)
+                                          reference_velocity=velocity_function2d)
 
 # ==========================================================================================
 # ==== Boundary

examples/fluid/pipe_flow_3d.jl

@@ -0,0 +1,45 @@
+# 3D channel flow simulation with open boundaries.
+
+using TrixiParticles
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.05
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 3
+
+# Make sure that the kernel support of fluid particles at an open boundary is always
+# 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 = 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
+
+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])
+
+flow_direction = [1.0, 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,
+              flow_direction=flow_direction, faces=(false, false, true, true, true, true),
+              tspan=tspan, smoothing_kernel=WendlandC2Kernel{3}(),
+              reference_velocity=velocity_function3d,
+              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],
+                         [domain_size[1], 0.0, domain_size[3]]))

examples/fsi/falling_spheres_2d.jl

@@ -63,7 +63,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
 
 # ==========================================================================================
 # ==== Boundary
-boundary_density_calculator = AdamiPressureExtrapolation()
+boundary_density_calculator = BernoulliPressureExtrapolation()
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              boundary_density_calculator,