From ef67faae5cf194e7af2da7879773793940375887 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 9 Mar 2026 16:47:21 +0100 Subject: [PATCH 1/2] correct language and general update pass --- README.md | 35 ++++++++++++++----------- docs/src/development.md | 14 +++++----- docs/src/examples.md | 4 +-- docs/src/general/neighborhood_search.md | 9 ++++--- docs/src/getting_started.md | 22 ++++++++-------- docs/src/gpu.md | 19 +++++++------- docs/src/index.md | 30 +++++++++++---------- docs/src/install.md | 23 ++++++++-------- docs/src/overview.md | 31 +++++++++++----------- docs/src/visualization.md | 16 ++++++----- 10 files changed, 106 insertions(+), 97 deletions(-) diff --git a/README.md b/README.md index 185f9317a0..27a54e37e6 100644 --- a/README.md +++ b/README.md @@ -15,30 +15,29 @@ TrixiP_logo

-**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** is a high-performance simulation framework for particle-based methods for complex multiphysics applications, written in [Julia](https://julialang.org). TrixiParticles.jl focuses on the following use cases: -- Accurate and efficient physics-based modelling of complex multiphysics problems. +- Accurate and efficient physics-based modeling 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. +- Accessible simulation setup 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). +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) ## Features -- Incompressible Navier-Stokes - - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), - Entropically Damped Artificial Compressibility (EDAC), - Implicit Incompressible SPH (IISPH) +- Incompressible Navier-Stokes flows + - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC), + Implicit Incompressible Smoothed Particle Hydrodynamics (IISPH) - Models: Surface Tension, Open Boundaries -- Solid-body mechanics - - Methods: Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM) +- Structural mechanics + - Methods: Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM) - Fluid-Structure Interaction -- Particle sampling of complex geometries from `.stl` and `.asc` files. +- Particle sampling of complex geometries from `.stl`, `.asc`, and `.dxf` files. - Output formats: - VTK -- Support for GPUs by Nvidia, AMD and Apple (experimental) +- GPU support for NVIDIA, AMD, and Apple devices ## Examples We provide several example simulation setups in the `examples` folder (which can be accessed from Julia via `examples_dir()`). @@ -68,7 +67,7 @@ We provide several example simulation setups in the `examples` folder (which can ## Installation -If you have not yet installed Julia, please [follow the instructions for your +If you have not installed Julia yet, please [follow the instructions for your operating system](https://julialang.org/downloads/platform/). TrixiParticles.jl works with Julia v1.10 and newer. We recommend using the latest stable release of Julia. @@ -115,7 +114,12 @@ Then start the simulation by executing julia> trixi_include(joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl")) ``` -This will open a new window with a 2D visualization of the final solution: +To visualize the result quickly, use Plots.jl: +```julia +julia> using Plots; plot(sol) +``` + +This opens a new window with a 2D visualization of the final solution: Further details can be found in the [documentation](https://trixi-framework.github.io/TrixiParticles.jl/stable). @@ -160,7 +164,7 @@ and ## Authors Erik Faulhaber (University of Cologne) and Niklas Neher (HLRS) implemented the foundations -for TrixiParticles.jl and are principal developers along with Sven Berger (hereon). +for TrixiParticles.jl and are principal developers along with Sven Berger (Hereon). The project was started by Michael Schlottke-Lakemper (University of Augsburg) and Gregor Gassner (University of Cologne), who provide scientific direction and technical advice. The full list of contributors can be found in [AUTHORS.md](AUTHORS.md). @@ -183,4 +187,3 @@ or [create an issue](https://github.com/trixi-framework/TrixiParticles.jl/issues

The project has benefited from funding from [hereon](https://www.hereon.de/), [HiRSE](https://www.helmholtz-hirse.de/), and through [ScienceServe](https://www.helmholtz.de/en/research/current-calls-for-applications/article/scienceserve-boosting-research-software-at-helmholtz/) for the MATRIX project. - diff --git a/docs/src/development.md b/docs/src/development.md index 8ffb75af52..09201863c1 100644 --- a/docs/src/development.md +++ b/docs/src/development.md @@ -3,24 +3,25 @@ ## Preview of the documentation -To generate the Documentation, first instantiate the `docs` environment -by executing the following command from the TrixiParticles.jl root directory: +To build the documentation, first instantiate the `docs` environment by running the +following command from the TrixiParticles.jl root directory: ```bash julia --project=docs -e "using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()" ``` -This command only has to be run once. After that, maintain the `docs` environment +You only need to run this command once. After that, maintain the `docs` environment as described under [Installation](@ref installation-issues). -With an instantiated `docs` environment, generate the docs with the following command (again from the TrixiParticles.jl root directory): +Once the `docs` environment is instantiated, build the documentation with the following +command (again from the TrixiParticles.jl root directory): ```bash julia --project=docs --color=yes docs/make.jl ``` -You can then open the generated files in `docs/build` with your webbrowser. +You can then open the generated files in `docs/build` in your web browser. Alternatively, run ```bash python3 -m http.server -d docs/build ``` -and open `localhost:8000` in your webbrowser. +and open `localhost:8000` in your web browser. ## Release management @@ -65,4 +66,3 @@ To create a new release for TrixiParticles.jl, perform the following steps: `-dev` suffix added. For example, if you just released `v0.3.0`, the new development version should be `v0.3.1-dev`. If you just released `v0.2.4`, the new development version should be `v0.2.5-dev`. - diff --git a/docs/src/examples.md b/docs/src/examples.md index f19e623e9e..59e2bda55c 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -42,7 +42,7 @@ ``` -## Fluid Structure Interaction +## Fluid-Structure Interaction ### Dam Break with Elastic Plate (`fsi/dam_break_plate_2d.jl`) ```@raw html @@ -59,7 +59,7 @@ ``` -## Structure Mechanics +## Structural Mechanics ### Oscillating Beam (`solid/oscillating_beam_2d.jl`) ```@raw html diff --git a/docs/src/general/neighborhood_search.md b/docs/src/general/neighborhood_search.md index b639e4a05f..93019eaffc 100644 --- a/docs/src/general/neighborhood_search.md +++ b/docs/src/general/neighborhood_search.md @@ -1,13 +1,14 @@ # Neighborhood Search -The neighborhood search is the most essential component for performance. +The neighborhood search is one of the most performance-critical components. We provide several implementations in the package [PointNeighbors.jl](https://github.com/trixi-framework/PointNeighbors.jl). -See the docs of this package for an overview and a comparison of different implementations. +See the PointNeighbors.jl documentation for an overview and a comparison of the +different implementations. !!! note "Usage" - To run a simulation with a neighborhood search implementation, pass a template of the - neighborhood search to the constructor of the [`Semidiscretization`](@ref). + To run a simulation with a neighborhood search implementation, pass a neighborhood + search template to the constructor of the [`Semidiscretization`](@ref). A template is just an empty neighborhood search with search radius `0.0`. See [`copy_neighborhood_search`](@ref) and the examples below for more details. ```jldoctest semi_example; output=false, setup = :(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); system1 = fluid_system; system2 = boundary_system) diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index a68bb620c0..08ebfeb066 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -1,14 +1,14 @@ # [Getting started](@id getting_started) -If you have not installed TrixiParticles.jl, please follow the instructions given [here](install.md). +If you have not installed TrixiParticles.jl yet, please follow the instructions in [Installation](@ref installation). -In the following sections, we will give a short introduction. For a more thorough discussion, take a look at our [Tutorials](tutorial.md). +This page provides a short introduction. For a broader introduction, take a look at our [Tutorials](tutorial.md). ## Running an Example -The easiest way to run a simulation is to run one of our predefined example files. +The easiest way to start is to run one of the predefined example files. We will run the file `examples/fluid/hydrostatic_water_column_2d.jl`, which simulates a fluid resting in a rectangular tank. Since TrixiParticles.jl uses multithreading, you should start Julia with the flag `--threads auto` (or, e.g. `--threads 4` for 4 threads). -In the Julia REPL, first load the package TrixiParticles.jl. +In the Julia REPL, first load `TrixiParticles.jl`. ```jldoctest getting_started julia> using TrixiParticles ``` @@ -29,8 +29,8 @@ This will open a new window with a 2D visualization of the final solution: For more information about visualization, see [Visualization](visualization.md). ## Running other Examples -You can find a list of our other predefined examples under [Examples](examples.md). -Execute them as follows from the Julia REPL by replacing `subfolder` and `example_name` +You can find more predefined examples under [Examples](examples.md). +Run them from the Julia REPL by replacing `subfolder` and `example_name`: ```julia julia> trixi_include(joinpath(examples_dir(), "subfolder", "example_name.jl")) ``` @@ -38,18 +38,18 @@ julia> trixi_include(joinpath(examples_dir(), "subfolder", "example_name.jl")) ## Modifying an example You can pass keyword arguments to the function `trixi_include` to overwrite assignments in the file. -With `trixi_include`, we can overwrite variables defined in the example file to run a different simulation without modifying the example file. +With `trixi_include`, we can overwrite variables defined in the example file to run a different simulation without modifying the file itself. ```jldoctest getting_started; filter = r".*"s julia> trixi_include(joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), initial_fluid_size=(1.0, 0.5)) ``` -This for example, will change the fluid size from ``(0.9, 1.0)`` to ``(1.0, 0.5)``. +This, for example, changes the fluid size from ``(0.9, 1.0)`` to ``(1.0, 0.5)``. -To understand why, take a look into the file `hydrostatic_water_column_2d.jl` in the subfolder `fluid` inside the examples directory, which is the file that we executed earlier. +To understand why, take a look at the file `hydrostatic_water_column_2d.jl` in the `fluid` subdirectory of the examples directory, which is the file we executed earlier. You can see that the initial size of the fluid is defined in the variable `initial_fluid_size`, which we could overwrite with the `trixi_include` call above. Another variable that is worth experimenting with is `fluid_particle_spacing`, which controls the resolution of the simulation in this case. A lower value will increase the resolution and the runtime. -## Set up you first simulation from scratch +## Set Up Your First Simulation from Scratch See [Set up your first simulation](tutorials/tut_setup.md). -Find an overview over the available tutorials under [Tutorials](tutorial.md). +An overview of the available tutorials is available under [Tutorials](tutorial.md). diff --git a/docs/src/gpu.md b/docs/src/gpu.md index bfcebe24e0..0b14c12c70 100644 --- a/docs/src/gpu.md +++ b/docs/src/gpu.md @@ -1,9 +1,8 @@ # [GPU Support](@id gpu_support) -GPU support is still an experimental feature that is actively being worked on. Currently, the [`WeaklyCompressibleSPHSystem`](@ref), [`TotalLagrangianSPHSystem`](@ref) and [`WallBoundarySystem`](@ref) support GPU execution. -We have tested GPU support on Nvidia, AMD and Apple GPUs. +We have tested GPU support on NVIDIA, AMD, and Apple GPUs. Note that most Apple GPUs do not support `Float64`. See [below on how to run single precision simulations](@ref single_precision). @@ -42,14 +41,14 @@ semi = Semidiscretization(fluid_system, boundary_system, └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` -At this point, we should run the simulation and make sure that it still works and that -the bounding box is large enough. +At this point, 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 parallelization backend that is used for the simulation. -On an Nvidia GPU, we specify: +Then, we only need to specify the parallelization backend used for the simulation. +On an NVIDIA GPU, we specify: ```julia using CUDA semi = Semidiscretization(fluid_system, boundary_system, @@ -65,7 +64,7 @@ semi = Semidiscretization(fluid_system, boundary_system, ``` Now, 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. +and their neighbors are 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). ## Run an existing example file on the GPU @@ -84,7 +83,7 @@ Note that in `examples/fluid/dam_break_2d.jl`, we explicitly set so that we can use `trixi_include` to replace this value. To run this simulation on a GPU, simply update `parallelization_backend` to the backend -of the installed GPU. We can run this simulation on an Nvidia GPU as follows. +of the installed GPU. We can run this simulation on an NVIDIA GPU as follows. ```julia using CUDA trixi_include(joinpath(examples_dir(), "fluid", "dam_break_2d_gpu.jl"), parallelization_backend=CUDABackend()) @@ -94,7 +93,7 @@ For AMD GPUs, use using AMDGPU trixi_include(joinpath(examples_dir(), "fluid", "dam_break_2d_gpu.jl"), parallelization_backend=ROCBackend()) ``` -For Apple GPUs (which don't support double precision, see below), use +For Apple GPUs (which do not support double precision, see below), use ```julia using Metal trixi_include_changeprecision(Float32, @@ -110,7 +109,7 @@ which is significantly faster on most GPUs and required for many Apple GPUs. To run a simulation with single precision, all `Float64` literals in an example file must be converted to `Float32` (e.g. `0.0` to `0.0f0`). -TrixiParticles provides a function to automate this conversion: +TrixiParticles.jl provides a function to automate this conversion: ```@docs trixi_include_changeprecision ``` diff --git a/docs/src/index.md b/docs/src/index.md index 27c1e9d474..19b4a32b87 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,24 +1,26 @@ # TrixiParticles.jl -**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** is a high-performance simulation framework for particle-based methods in complex multiphysics applications. It combines an accessible user interface with an extensible architecture for developing new methods, while also providing GPU-accelerated 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. +- Development of new particle-based methods and models through an extensible architecture that is not tied to a single numerical method. +- Accurate, reliable, and efficient physics-based modeling of complex multiphysics problems through a flexible configuration system, high performance, and a broad set of validation and test cases. +- Accessible simulation setup for educational purposes, including student projects, coursework, and thesis work, supported by extensive documentation and readable configuration files. -Its features include: +Its main features include: ## Features -- Incompressible Navier-Stokes - - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC) - - Models: Surface Tension -- Solid-body mechanics - - Methods: Total Lagrangian SPH (TLSPH) +- Incompressible Navier-Stokes flows + - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC), Implicit Incompressible Smoothed Particle Hydrodynamics (IISPH) + - Models: Surface Tension, Open Boundaries +- Structural mechanics + - Methods: Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM) - Fluid-Structure Interaction +- Particle sampling of complex geometries from `.stl`, `.asc`, and `.dxf` files - Output formats: - VTK +- Experimental GPU support for NVIDIA, AMD, and Apple devices ## Examples ```@raw html @@ -42,14 +44,14 @@ Its features include: ``` -## Quickstart +## Quick Start 1. [Installation](@ref installation) 2. [Getting started](@ref getting_started) -If you have any questions concerning **TrixiParticles.jl** you can join our community [on Slack](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g) or open an issue with your question. +If you have questions about **TrixiParticles.jl**, join our community [on Slack](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g) or open an issue. -## Start with development -To get started with development have a look at these pages: +## Getting Started with Development +If you want to contribute or extend the code, start with: 1. [Installation](@ref installation) 2. [Development](@ref development) diff --git a/docs/src/install.md b/docs/src/install.md index 15c6668b9c..ad379e6859 100644 --- a/docs/src/install.md +++ b/docs/src/install.md @@ -1,7 +1,7 @@ # [Installation](@id installation) ## Setting up Julia -If you have not yet installed Julia, please [follow the instructions on the +If you have not installed Julia yet, please [follow the instructions on the official website](https://julialang.org/downloads/). TrixiParticles.jl works with Julia v1.10 and newer. We recommend using the latest stable release of Julia. @@ -35,26 +35,25 @@ julia --project=run ``` from the TrixiParticles.jl root directory. -The advantage of using a separate `run` directory is that you can also add other -related packages (e.g., OrdinaryDiffEq.jl, see above) to the project in the `run` folder -and always have a reproducible environment at hand to share with others. +The advantage of using a separate `run` directory is that you can add other related +packages (e.g., OrdinaryDiffEq.jl, see above) to the project in that folder while +keeping a reproducible environment that is easy to share with others. ## Optional software/packages -- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- A Julia package of ordinary differential equation solvers that is used in the examples -- [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia Plotting library that is used in some examples +- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- Julia package of ordinary differential equation solvers used in the examples +- [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia plotting library used in some examples - [PythonPlot.jl](https://github.com/JuliaPy/PythonPlot.jl) -- Plotting library that can be used instead of Plots.jl -- [ParaView](https://www.paraview.org/) -- Software that can be used for visualization of results +- [ParaView](https://www.paraview.org/) -- Visualization software for simulation results ## [Common issues](@id installation-issues) -If you followed the [installation instructions for developers](@ref for-developers) and you -run into any problems with packages when pulling the latest version of TrixiParticles.jl, -start Julia with the project in the `run` folder, +If you followed the [installation instructions for developers](@ref for-developers) and run +into package issues after pulling the latest version of TrixiParticles.jl, start Julia with +the project in the `run` folder, ```bash julia --project=run ``` -update all packages in that project, resolve all conflicts in the project, and install all -new dependencies: +then update packages, resolve dependency conflicts, and install new dependencies: ```julia julia> using Pkg diff --git a/docs/src/overview.md b/docs/src/overview.md index a1b6a70d5a..fc500ab85d 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -1,18 +1,18 @@ # 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 API reference is not collected on a single page. +Instead, it is split into sections that largely mirror the source tree. +These sections combine API documentation with short explanations of the +underlying 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 -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. +To run a simulation, TrixiParticles.jl solves an ordinary differential equation, +typically with a time integration scheme from OrdinaryDiffEq.jl. These schemes are +used to integrate ``\mathrm{d}u/\mathrm{d}t`` and ``\mathrm{d}v/\mathrm{d}t``, +where ``u`` represents particle positions and ``v`` particle properties such as +velocity and density. During a single time step or an intermediate step of the time integration scheme, the functions `drift!` and `kick!` are invoked, followed by the functions depicted in this diagram (with key parts highlighted in orange/yellow). @@ -78,9 +78,10 @@ flowchart TD ``` ## Structure -What we refer to as schemes are various models such as Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH) -or Total Lagrangian Smoothed Particle Hydrodynamics (TLSPH). These schemes are categorized based on the applicable -physical regimes, namely fluid, solid, gas, and others. Each scheme comprises at least two files: a `system.jl` file -and an `rhs.jl` file. The `system.jl` file provides the data structure holding the particles of this scheme and some -routines, particularly those for allocation and the main update routines, excluding system interactions. -The interactions between particles of this scheme (and with particles of other schemes) are handled in the `rhs.jl` file. +In the codebase, a scheme denotes a particle method or model such as Weakly Compressible +SPH (WCSPH) or Total Lagrangian SPH (TLSPH). Schemes are organized by application area, +for example fluid, structure, and boundary systems. A scheme typically comprises at least +two files: a `system.jl` file and an `rhs.jl` file. The `system.jl` file defines the data +structure that stores the particles of the scheme together with routines for allocation and +main updates that do not involve particle interactions. The `rhs.jl` file contains the +interaction terms between particles of the same scheme and between different schemes. diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 1f74108b4f..9dda94bfba 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -2,8 +2,10 @@ ## Export VTK files You can export particle data as VTK files by using the [`SolutionSavingCallback`](@ref). -All our [predefined examples](examples.md) are already using this callback to export VTK files to the `out` directory (relative to the directory that you are running Julia from). -VTK files can be read by visualization tools like [ParaView](https://www.paraview.org/) and [VisIt](https://visit.llnl.gov/). +All [predefined examples](examples.md) already use this callback to export VTK files to the `out` +directory relative to the current working directory. +VTK files can be opened in visualization tools such as [ParaView](https://www.paraview.org/) +and [VisIt](https://visit.llnl.gov/). ### ParaView @@ -29,8 +31,8 @@ Then, in the Properties panel (bottom left), adjust the following settings: ![image](https://github.com/user-attachments/assets/194d9a09-5937-4ee4-b229-07078afe3ff0) #### Visualization with Macro -To simplify the visualization of your particle data in ParaView, you can use a macro. -This macro automates the manual steps in the previous section to a single click of a button. +To simplify visualization of particle data in ParaView, you can use a macro. +It reduces the manual steps from the previous section to a single click. Install the macro as follows. 1. **Save the macro code** (see below) as a `.py` file, e.g. `PointGaussianMacro.py`. @@ -77,8 +79,10 @@ sourceDisplay.GaussianRadius = 0.5 ``` #### Show results -To now view the result variables **first** make sure you have "fluid_1.pvd" highlighted in the "Pipeline Browser" then select them in the variable selection combo box (see picture below). -Let's, for example, pick "density". To now view the time progression of the result hit the "play button" (see picture below). +To view the result variables, first make sure that "fluid_1.pvd" is highlighted in the +"Pipeline Browser", then select a variable in the variable-selection combo box +(see the image below). For example, choose "density". To view the time evolution, +press the play button (also shown below). ![image](https://github.com/user-attachments/assets/10dcf7eb-5808-4d4d-9db8-4beb25b5e51a) ## API From a84e50c792a2d99e872589054bc0235f906cef9d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 9 Mar 2026 18:43:51 +0100 Subject: [PATCH 2/2] improve maths --- docs/src/systems/boundary.md | 53 +++++--- docs/src/systems/entropically_damped_sph.md | 37 ++++-- docs/src/systems/fluid.md | 117 +++++++++++------ .../systems/implicit_incompressible_sph.md | 100 +++++++------- docs/src/systems/total_lagrangian_sph.md | 8 +- docs/src/systems/weakly_compressible_sph.md | 123 ++++++++++++++---- docs/src/time_integration.md | 4 +- 7 files changed, 291 insertions(+), 151 deletions(-) diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 94fbe6ce05..46463b10ac 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -16,6 +16,11 @@ Pages = [joinpath("schemes", "boundary", "prescribed_motion.jl")] # [Boundary Models](@id boundary_models) +!!! note + The pairwise interaction terms below are written in force form, following the SPH literature. + TrixiParticles.jl applies the corresponding accelerations internally. Where the implemented + boundary discretization differs from the literature formula, both forms are stated explicitly. + ## Dummy Particles Boundaries modeled as dummy particles, which are treated like fluid particles, @@ -25,7 +30,7 @@ dummy particles need to have a mass corresponding to the fluid's rest density, w "hydrodynamic mass", as opposed to mass corresponding to the material density of a [`TotalLagrangianSPHSystem`](@ref). -Here, `initial_density` and `hydrodynamic_mass` are vectors that contains the initial density +Here, `initial_density` and `hydrodynamic_mass` are vectors that contain the initial density and the hydrodynamic mass respectively for each boundary particle. Note that when used with [`SummationDensity`](@ref) (see below), this is only used to determine the element type and the number of boundary particles. @@ -37,14 +42,16 @@ This should be the same as for the adjacent fluid system with the largest smooth In the literature, this kind of boundary particles is referred to as "dummy particles" ([Adami et al., 2012](@cite Adami2012) and [Valizadeh & Monaghan, 2015](@cite Valizadeh2015)), -"frozen fluid particles" ([Akinci et al., 2012](@cite Akinci2012)) or "dynamic boundaries [Crespo et al., 2007](@cite Crespo2007). +"frozen fluid particles" ([Akinci et al., 2012](@cite Akinci2012)) or "dynamic boundaries" ([Crespo et al., 2007](@cite Crespo2007)). The key detail of this boundary condition and the only difference between the boundary models in these references is the way the density and pressure of boundary particles is computed. -Since boundary particles are treated like fluid particles, the force -on fluid particle ``a`` due to boundary particle ``b`` is given by +For the standard summation-density pressure force, the force on fluid particle ``a`` +due to boundary particle ``b`` is ```math -f_{ab} = m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_{r_a} W(\Vert r_a - r_b \Vert, h). +\bm{f}_{ab}^{p} += -m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h). ``` The quantities to be defined here are the density ``\rho_b`` and pressure ``p_b`` of the boundary particle ``b``. @@ -99,11 +106,16 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den ``` #### 2. [`BernoulliPressureExtrapolation`](@ref) -Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation: +Identical to the pressure ``p_b`` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds an additional dynamic pressure term. For moving wall boundaries, the implementation uses +```math +p_b = \frac{\sum_f (p_f + p_{f,\mathrm{dyn}} + \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)}, +``` +with ```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)} +p_{f,\mathrm{dyn}} = \frac{1}{2} \, \text{factor} \, \rho_f +\frac{\left((\bm{v}_b - \bm{v}_f) \cdot \bm{r}_{bf}\right)^2}{\Vert \bm{r}_{bf} \Vert}, ``` -where ``\mathbf{v}_f`` is the velocity of the fluid and ``\mathbf{v}_{\text{body}}`` is the velocity of the body. +where ``\bm{v}_f`` is the fluid velocity and ``\bm{v}_b`` is the boundary velocity. This adjustment provides a higher boundary pressure for bodies moving with a relative velocity to the fluid to prevent penetration. This modification is original and not derived from any literature source. @@ -123,17 +135,16 @@ reference pressure (the corresponding pressure to the reference density by the s #### 6. [`PressureMirroring`](@ref) Instead of calculating density and pressure for each boundary particle, we modify the -momentum equation, +pressure force, ```math -\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab} +\bm{F}_a^{p} = -m_a \sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab}, ``` to replace the unknown density $\rho_b$ if $b$ is a boundary particle by the reference density and the unknown pressure $p_b$ if $b$ is a boundary particle by the pressure $p_a$ of the -interacting fluid particle. -The momentum equation therefore becomes +interacting fluid particle. The force therefore becomes ```math -\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af} --\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab}, +\bm{F}_a^{p} = -m_a \sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af} +-m_a \sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab}, ``` where the first sum is over all fluid particles and the second over all boundary particles. @@ -168,18 +179,22 @@ condition is applied. ## Repulsive Particles -Boundaries modeled as boundary particles which exert forces on the fluid particles ([Monaghan, Kajtar, 2009](@cite Monaghan2009)). -The force on fluid particle ``a`` due to boundary particle ``b`` is given by +Boundaries modeled as boundary particles which exert repulsive interactions on the fluid particles ([Monaghan, Kajtar, 2009](@cite Monaghan2009)). +The literature force on fluid particle ``a`` due to boundary particle ``b`` is ```math -f_{ab} = m_a \left(\tilde{f}_{ab} - m_b \Pi_{ab} \nabla_{r_a} W(\Vert r_a - r_b \Vert, h)\right) +\bm{f}_{ab} = m_a \left(\tilde{\bm{f}}_{ab} - m_b \Pi_{ab} +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h)\right) ``` with ```math -\tilde{f}_{ab} = \frac{K}{\beta^{n-1}} \frac{r_{ab}}{\Vert r_{ab} \Vert (\Vert r_{ab} \Vert - d)} \Phi(\Vert r_{ab} \Vert, h) +\tilde{\bm{f}}_{ab} = +\frac{K}{\beta^{n-1}} \frac{\bm{r}_{ab}} +{\Vert \bm{r}_{ab} \Vert (\Vert \bm{r}_{ab} \Vert - d)} +\Phi(\Vert \bm{r}_{ab} \Vert, h)\, \frac{2 m_b}{m_a + m_b}, ``` where ``m_a`` and ``m_b`` are the masses of fluid particle ``a`` and boundary particle ``b`` -respectively, ``r_{ab} = r_a - r_b`` is the difference of the coordinates of particles +respectively, ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` is the difference of the coordinates of particles ``a`` and ``b``, ``d`` denotes the boundary particle spacing and ``n`` denotes the number of dimensions (see [Monaghan & Kajtar, 2009](@cite Monaghan2009), Equation (3.1) and [Valizadeh & Monaghan, 2015](@cite Valizadeh2015)). Note that the repulsive acceleration $\tilde{f}_{ab}$ does not depend on the masses of diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 96acbad352..be1772bcd7 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -3,42 +3,55 @@ As opposed to the [weakly compressible SPH scheme](weakly_compressible_sph.md), which uses an equation of state, this scheme uses a pressure evolution equation to calculate the pressure ```math -\frac{\mathrm{d} p_a}{\mathrm{d}t} = - \rho c_s^2 \nabla \cdot v + \nu \nabla^2 p, +\frac{\mathrm{d} p_a}{\mathrm{d}t} = - \rho_a c_s^2 (\nabla \cdot v)_a + \nu_{\mathrm{EDAC}} (\nabla^2 p)_a, ``` which is derived by [Clausen (2013)](@cite Clausen2013). This equation is similar to the continuity equation (first term, see -[`ContinuityDensity`](@ref)), but also contains a pressure damping term (second term, similar to density diffusion +[`ContinuityDensity`](@ref)), but also contains a pressure damping term (second term, similar to density diffusion, see [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion)), which reduces acoustic pressure waves through an entropy-generation mechanism. -The pressure evolution is discretized with the SPH method by [Ramachandran (2019)](@cite Ramachandran2019) as following: +The pressure evolution is discretized with the SPH method by [Ramachandran (2019)](@cite Ramachandran2019) as follows: The first term is equivalent to the classical artificial compressible methods, which are commonly motivated by assuming the artificial equation of state ([`StateEquationCole`](@ref) with `exponent=1`) and is discretized as ```math -- \rho c_s^2 \nabla \cdot v = \sum_{b} m_b \frac{\rho_a}{\rho_b} c_s^2 v_{ab} \cdot \nabla_{r_a} W(\Vert r_a - r_b \Vert, h), +\left.- \rho c_s^2 \nabla \cdot v \right|_a += \sum_{b} m_b \frac{\rho_a}{\rho_b} c_s^2 v_{ab} \cdot \nabla_{r_a} W(\Vert r_a - r_b \Vert, h), ``` where ``\rho_a``, ``\rho_b``, ``r_a``, ``r_b``, denote the density and coordinates of particles ``a`` and ``b`` respectively, ``c_s`` is the speed of sound and ``v_{ab} = v_a - v_b`` is the difference in the velocity. The second term smooths the pressure through the introduction of entropy and is discretized as ```math -\nu \nabla^2 p = \frac{V_a^2 + V_b^2}{m_a} \tilde{\eta}_{ab} \frac{p_{ab}}{\Vert r_{ab}^2 \Vert + \eta h_{ab}^2} \nabla_{r_a} -W(\Vert r_a - r_b \Vert, h) \cdot r_{ab}, +\left.\nu_{\mathrm{EDAC}} \nabla^2 p \right|_a += \sum_b \frac{V_a^2 + V_b^2}{m_a}\, +\tilde{\eta}_{ab}\, +\frac{p_{ab}}{\Vert r_{ab} \Vert^2 + 0.01 h_{ab}^2}\, +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h) \cdot r_{ab}, ``` -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively and ``p_{ab}= p_a -p_b`` is the difference in the pressure. +where ``V_a``, ``V_b`` denote the particle volumes, ``p_{ab}= p_a - p_b``, +``r_{ab} = r_a - r_b``, and ``h_{ab} = \frac{1}{2}(h_a + h_b)``. -The viscosity parameter ``\eta_a`` for a particle ``a`` is given as +The dynamic EDAC viscosity for particle ``a`` is ```math -\eta_a = \rho_a \frac{\alpha h c_s}{8}, +\eta_a = \rho_a \nu_{\mathrm{EDAC}}, ``` -where it is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5`` +with +```math +\nu_{\mathrm{EDAC}} = \frac{\alpha h c_s}{8}, +``` +and the harmonic mean +```math +\tilde{\eta}_{ab} = \frac{2 \eta_a \eta_b}{\eta_a + \eta_b}. +``` +It is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5`` is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). !!! note > The EDAC formulation keeps the density constant and this eliminates the need for the continuity equation - > or the use of a summation density to find the pressure. However, in SPH discretizations, ``m/\rho`` - > is typically used as a proxy for the particle volume. The density of the fluids can + > or the use of a summation density to find the pressure. However, in SPH discretizations, ``m/\rho`` + > is typically used as a proxy for the particle volume. The density of the fluids can > therefore be computed using the summation density approach. [Ramachandran2019](@cite) diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..a22bf11e9e 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -4,6 +4,12 @@ Currently available fluid methods are the [weakly compressible SPH method](@ref [entropically damped artificial compressibility for SPH](@ref edac). This page lists models and techniques that apply to both of these methods. +!!! note + The formulas on this page follow the force notation commonly used in the SPH literature. + TrixiParticles.jl usually evaluates the corresponding acceleration contributions internally. + Whenever the implemented discretization differs from the literature formula, both forms are + stated explicitly. + ## [Viscosity](@id viscosity_sph) Viscosity is a critical physical property governing momentum diffusion within a fluid. @@ -55,10 +61,10 @@ by Balsara ([Balsara1995](@cite)) or Morris ([Morris1997](@cite)). ##### Mathematical Formulation -The force exerted by particle ``b`` on particle ``a`` due to artificial viscosity is given by: +The force exerted by particle ``b`` on particle ``a`` due to artificial viscosity is given by ```math -F_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab} +\bm{F}_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla_a W_{ab}. ``` where: @@ -75,18 +81,19 @@ where: - ``c`` is the local speed of sound, - ``\bar{\rho}_{ab}`` is the arithmetic mean of the densities of particles ``a`` and ``b``. -The term ``\mu_{ab}`` is defined as: +The term ``\mu_{ab}`` is defined as ```math -\mu_{ab} = \frac{h \, v_{ab} \cdot r_{ab}}{\Vert r_{ab} \Vert^2 + \epsilon h^2}, +\mu_{ab} = \frac{h \, \bm{v}_{ab} \cdot \bm{r}_{ab}} + {\Vert \bm{r}_{ab} \Vert^2 + \epsilon h^2}, ``` with: - ``h`` being the smoothing length, - ``\epsilon`` a small parameter to prevent singularities, -- ``r_{ab} = r_a - r_b`` representing the difference of the coordinate vectors, -- ``v_{ab} = v_a - v_b`` representing the relative velocity between particles. +- ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` representing the difference of the coordinate vectors, +- ``\bm{v}_{ab} = \bm{v}_a - \bm{v}_b`` representing the relative velocity between particles. ##### Resolution Dependency and Effective Viscosity @@ -109,19 +116,21 @@ This results in a more realistic representation of flow dynamics in weakly compr ##### Mathematical Formulation -An additional force term ``\tilde{f}_{ab}`` is introduced to the pressure gradient force ``f_{ab}`` between particles ``a`` and ``b``: +An additional force term ``\tilde{\bm{F}}_{ab}`` is introduced in the momentum equation: ```math -\tilde{f}_{ab} = m_a m_b \frac{(\mu_a + \mu_b)\, r_{ab} \cdot \nabla W_{ab}}{\rho_a \rho_b (\Vert r_{ab} \Vert^2 + \epsilon h^2)}\, v_{ab}, +\tilde{\bm{F}}_{ab} = +m_a m_b \frac{(\mu_a + \mu_b)\, \bm{r}_{ab} \cdot \nabla_a W_{ab}} +{\rho_a \rho_b (\Vert \bm{r}_{ab} \Vert^2 + \epsilon h^2)}\, \bm{v}_{ab}, ``` where: -- ``\mu_a = \rho_a \nu`` and ``\mu_b = \rho_b \nu`` represent the dynamic viscosities of particles ``a``and ``b`` (with ``\nu`` being the kinematic viscosity), -- ``r_{ab} = r_a - r_b`` represents the difference of the coordinate vectors, -- ``v_{ab} = v_a - v_b`` represents the relative velocity between particles. +- ``\mu_a = \rho_a \nu`` and ``\mu_b = \rho_b \nu`` represent the dynamic viscosities of particles ``a`` and ``b`` (with ``\nu`` being the kinematic viscosity), +- ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` represents the difference of the coordinate vectors, +- ``\bm{v}_{ab} = \bm{v}_a - \bm{v}_b`` represents the relative velocity between particles, - `` h `` is the smoothing length, -- `` \nabla W_{ab} `` is the gradient of the smoothing kernel, +- `` \nabla_a W_{ab} `` is the gradient of the smoothing kernel, - `` \epsilon `` is a small parameter to prevent singularities. #### ViscosityAdami @@ -132,19 +141,24 @@ while minimizing compressibility effects. This results in accurate laminar flow ##### Mathematical Formulation -The viscous interaction is modeled through a shear force for incompressible flows: +The viscous interaction is modeled through the following pairwise force: ```math -f_{ab} = \sum_w \bar{\eta}_{ab} \left( V_a^2 + V_b^2 \right) \frac{v_{ab}}{||r_{ab}||^2 + \epsilon h_{ab}^2} \, (\nabla W_{ab} \cdot r_{ab}), +\bm{F}_{ab}^{\nu} = +\left( V_a^2 + V_b^2 \right)\, +\bar{\eta}_{ab}\, +\frac{\nabla_a W_{ab} \cdot \bm{r}_{ab}} +{\Vert \bm{r}_{ab} \Vert^2 + \epsilon h_{ab}^2}\, +\bm{v}_{ab}. ``` where: -- `` r_{ab} = r_a - r_b `` is the difference of the coordinate vectors, -- `` v_{ab} = v_a - v_b `` is their relative velocity, +- `` \bm{r}_{ab} = \bm{r}_a - \bm{r}_b `` is the difference of the coordinate vectors, +- `` \bm{v}_{ab} = \bm{v}_a - \bm{v}_b `` is their relative velocity, - `` V_a = m_a / \rho_a`` and `` V_b = m_b / \rho_b`` are the particle volumes, -- `` h_{ab} `` is the smoothing length, -- `` \nabla W_{ab} `` is the gradient of the smoothing kernel, +- `` h_{ab} = \frac{1}{2}(h_a + h_b) `` is the arithmetic mean of the smoothing lengths, +- `` \nabla_a W_{ab} `` is the gradient of the smoothing kernel, - `` \epsilon `` is a small parameter that prevents singularities (see [Ramachandran (2019)](@cite Ramachandran2019)). The inter-particle-averaged shear stress is defined as: @@ -225,10 +239,10 @@ The surface normal at a particle is derived from the color field, a scalar field to distinguish between different fluid phases or between fluid and air. The color field gradients point towards the interface, and the normalized gradient defines the surface normal direction. -The simplest SPH formulation for a surface normal, ``n_a`` is given as +In the literature, the unnormalized surface normal ``\bm{n}_a`` is commonly written as ```math -n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, +\bm{n}_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, ``` where: @@ -238,12 +252,19 @@ where: - ``\rho_b`` is the density of particle ``b``, - ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``. +Implementation note: for the single-fluid surface-normal calculation in TrixiParticles.jl, +this reduces to +```math +\bm{n}_a = \sum_b m_b \frac{1}{\rho_b} \nabla_a W_{ab}, +``` +i.e. effectively ``c_b = 1`` for neighboring fluid particles. + #### Normalization of surface normals The calculated normals are normalized to unit vectors: ```math -\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}. +\hat{\bm{n}}_a = \frac{\bm{n}_a}{\Vert \bm{n}_a \Vert}. ``` Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. @@ -303,10 +324,10 @@ It is defined by the distance between particles and the support radius ``h_c``, - Particles within half the support radius experience a repulsive force to prevent clustering. - Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion. -Mathematically: +The pairwise cohesion force is ```math -F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, +\bm{F}_{\text{cohesion}} = -\sigma m_b C(r) \frac{\bm{r}}{\Vert \bm{r} \Vert}, ``` where ``C(r)``, the cohesion kernel, is defined as: @@ -315,21 +336,30 @@ where ``C(r)``, the cohesion kernel, is defined as: C(r)=\frac{32}{\pi h_c^9} \begin{cases} (h_c-r)^3 r^3, & \text{if } 2r > h_c, \\ -2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ +2(h_c-r)^3 r^3 - \frac{h_c^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ 0, & \text{otherwise.} \end{cases} ``` #### Surface area minimization force -The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. -It acts based on the difference in surface normals: +The surface area minimization force models curvature reduction and acts on the +difference in surface normals: ```math -F_{\text{curvature}} = -\sigma (n_a - n_b), +\bm{F}_{\text{curvature}} = -\sigma (\bm{n}_a - \bm{n}_b), ``` -where ``n_a`` and ``n_b`` are the surface normals of the interacting particles. +where ``\bm{n}_a`` and ``\bm{n}_b`` are the surface normals of the interacting particles. +Implementation note: TrixiParticles.jl uses +```math +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d} t}\right|_{ab}^{\text{curvature}} += -\sigma h (\bm{n}_a - \bm{n}_b), +``` +which uses the same normal-difference direction but scales the magnitude with the +local smoothing length ``h``. For constant ``h``, this factor can be absorbed into +the coefficient ``\sigma``; for variable smoothing lengths, it makes the curvature +contribution depend on the local kernel support. #### Wall adhesion force @@ -337,7 +367,7 @@ This force models the interaction between fluid and solid boundaries, simulating It uses a custom kernel with a peak at 0.75 times the support radius: ```math -F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, +\bm{F}_{\text{adhesion}} = -\beta m_b A(r) \frac{\bm{r}}{\Vert \bm{r} \Vert}, ``` where ``A(r)`` is the adhesion kernel: @@ -358,11 +388,13 @@ The method described by [Morris](@cite Morris2000) estimates curvature by combin The computed curvature is then used to determine forces acting perpendicular to the interface. While this method provides accurate surface tension forces, it does not explicitly conserve momentum. -In the Morris model, surface tension is computed based on local interface curvature ``\kappa`` and the unit surface normal ``\hat{n}.`` -By estimating ``\hat{n}`` and ``\kappa`` at each particle near the interface, the surface tension force for particle a can be written as: +In the Morris model, surface tension is computed based on local interface curvature ``\kappa`` and the unit surface normal ``\hat{\bm{n}}``. +By estimating ``\hat{\bm{n}}`` and ``\kappa`` at each particle near the interface, the +surface tension force for particle ``a`` can be written as ```math -F_{\text{surface tension}} = - \sigma \frac{\kappa_a}{\rho_a}\hat{n}_a +\bm{F}_{a}^{\sigma} += - m_a \sigma \frac{\kappa_a}{\rho_a}\hat{\bm{n}}_a. ``` This formulation focuses directly on geometric properties of the interface, making it relatively straightforward to implement when a reliable interface detection @@ -379,29 +411,38 @@ where accumulated numerical error can be significant. #### Stress tensor formulation -The surface tension force can be seen as a divergence of a stress tensor ``S`` +The surface tension force can be written as the divergence of a stress tensor ``\bm{S}``: ```math -F_{\text{surface tension}} = \nabla \cdot S, +\bm{F}_{a}^{\sigma} = m_a \nabla \cdot \bm{S}, ``` -with ``S`` defined as +with ```math -S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}), +\bm{S} = \sigma \delta_s (I - \hat{\bm{n}} \otimes \hat{\bm{n}}). ``` with: - ``\delta_s``: Surface delta function, -- ``\hat{n}``: Unit normal vector, +- ``\hat{\bm{n}}``: Unit normal vector, - ``I``: Identity matrix. 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} +\bm{F}_{a}^{\sigma} += m_a \sum_b \frac{m_b}{\rho_a \rho_b} (\bm{S}_a + \bm{S}_b) \nabla_a W_{ab}. +``` + +Implementation note: TrixiParticles.jl evaluates the corresponding acceleration and uses +the stabilized stress tensor +```math +\bm{S}_a^{\text{impl}} += \delta_{s,a} (I - \hat{\bm{n}}_a \otimes \hat{\bm{n}}_a) - \delta_{s,\max} I, ``` +with the factor ``\sigma`` applied outside the pairwise sum. #### Advantages and limitations diff --git a/docs/src/systems/implicit_incompressible_sph.md b/docs/src/systems/implicit_incompressible_sph.md index 0927aa8821..575d7d84a9 100644 --- a/docs/src/systems/implicit_incompressible_sph.md +++ b/docs/src/systems/implicit_incompressible_sph.md @@ -5,7 +5,10 @@ is a method that achieves incompressibility by solving the pressure Poisson equa The resulting linear system is iteratively solved with the relaxed Jacobi method. Unlike the [weakly compressible SPH method](@ref wcsph), incompressible methods determine pressure by enforcing the incompressibility constraint rather than using an equation of -state. +state. In the derivation below, we keep the force notation of +[Ihmsen et al. (2013)](@cite Ihmsen2013), which is also the standard presentation in the +IISPH literature. Internally, TrixiParticles.jl applies the equivalent pressure +accelerations ``\bm{F}_i^p / m_i``; this changes only the presentation, not the algebra. ```@autodocs Modules = [TrixiParticles] @@ -28,14 +31,15 @@ difference yields The divergence in the right-hand side is discretized with the SPH discretization for particle ``i`` as ```math --\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \nabla W_{ij}, +-\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \cdot \nabla W_{ij}, ``` where ``\bm{v}_{ij} = \bm{v}_i - \bm{v}_j``. Together, the following discretized version of the continuity equation for a particle ``i`` is achieved: ```math -\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij}. +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} += \sum_j m_j \bm{v}_{ij}(t+\Delta t) \cdot \nabla W_{ij}. ``` Note that the linear system is only solved for fluid particles, so ``i`` always represents @@ -50,12 +54,12 @@ Using the semi-implicit Euler method, we can obtain the velocity in the next tim ``` where ``\bm{F}_i^{\text{adv}}`` denotes all non-pressure forces such as gravity, viscosity, surface -tension and more, while ``\bm{F}_i^p``denotes the unknown pressure forces, which we +tension and more, while ``\bm{F}_i^p`` denotes the unknown pressure forces, which we want to solve for. Note that the IISPH is an incompressible method, which means that the density of the -fluid remain constant over time. By assuming a fixed reference density ``\rho_0`` for all -fluid particle over the whole time of the simulation, the density value at the next time +fluid remains constant over time. By assuming a fixed reference density ``\rho_0`` for all +fluid particles over the whole simulation, the density value at the next time step ``\rho_i(t + \Delta t)`` also has to be this rest density. So ``\rho_0`` can be plugged in for ``\rho_i(t + \Delta t)`` in the equation above. @@ -72,7 +76,7 @@ Using this predicted velocity and the continuity equation, a predicted density c in a similar way as ```math -\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t). +\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}}(t+\Delta t) \cdot \nabla W_{ij}(t). ``` To achieve the rest density, the unknown pressure forces must counteract the compression @@ -81,7 +85,7 @@ the predicted density and the reference density. Therefore, the following equation needs to be fulfilled: ```math -\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. +\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \cdot \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. ``` This expression is derived by substituting the reference density ``\rho_0`` for @@ -130,10 +134,10 @@ The pressure acceleration is given by: The ``d_{ii}p_i`` value describes the displacement of particle ``i`` because of the particle ``i`` and ``d_{ij}p_j`` describes the influence from the neighboring particles ``j``. -Using this new values the linear system can be rewritten as +Using these values, the linear system can be rewritten as ```math -\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \nabla W_{ij}, +\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \cdot \nabla W_{ij}, ``` where the first sum over ``k`` loops over all neighbor particles of ``i`` and @@ -150,21 +154,21 @@ To separate this sum, it can be written as With this separation, the equation for the linear system can again be rewritten as ```math -\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \nabla W_{ij}. +\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji}) \cdot \nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \cdot \nabla W_{ij}. ``` In this formulation all coefficients that are getting multiplied with the pressure value ``p_i`` are separated from the other. The diagonal elements ``a_{ii}`` can therefore be defined as: ```math -a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij}. +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji}) \cdot \nabla W_{ij}. ``` The remaining part of the equation represents the influence of the other pressure values ``p_j``. ​Hence, the final relaxed Jacobi iteration takes the form: ```math -p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right). +p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \cdot \nabla W_{ij} \right). ``` Because interactions are local, limited to particles within the kernel support defined by @@ -209,7 +213,7 @@ as only isolated or almost isolated particles are affected. ## Boundary Handling The previously introduced formulation did not distinguish between fluid and boundary -particles. To account boundary interactions correctly, a few modifications to the previous +particles. To account for boundary interactions correctly, a few modifications to the previous equations are required. First, the discretized form of the continuity equation must be adapted for the case in which @@ -220,7 +224,7 @@ neighboring fluid particles (indexed by ``f``) and neighboring boundary particle The updated discretized continuity equation becomes: ```math -\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_f m_f \bm{v}_{if}(t+\Delta t) \nabla W_{if} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \nabla W_{ib}. +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_f m_f \bm{v}_{if}(t+\Delta t) \cdot \nabla W_{if} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \cdot \nabla W_{ib}. ``` Since boundary particles have zero velocity, the difference between the fluid @@ -229,13 +233,13 @@ particle's velocity ``\bm{v}_{ib}(t+\Delta t) = \bm{v}_{i}(t+\Delta t)``. Accordingly, the predicted density ``\rho^{\text{adv}}`` becomes: ```math -\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_f m_f \bm{v}_{if}^{\text{adv}} \nabla W_{if}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \nabla W_{ib}(t). +\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_f m_f \bm{v}_{if}^{\text{adv}} \cdot \nabla W_{if}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \cdot \nabla W_{ib}(t). ``` This leads to the following updated formulation of the linear system: ```math -\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \nabla W_{if} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \cdot \nabla W_{if} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \cdot \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. ``` Note that, since boundary particles are fixed, the force ``F_b^p`` is zero and does not appear @@ -244,12 +248,12 @@ in this equation. The pressure force acting on a fluid particle is computed as: ```math -\bm{F}_i^p(t) = -\sum_f m_f \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_f(t)}{\rho_f^2(t)} \right) \nabla W_{if}(t) - \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). +\bm{F}_i^p(t) = -m_i \sum_f m_f \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_f(t)}{\rho_f^2(t)} \right) \nabla W_{if}(t) - m_i \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). ``` This also leads to an updated version of the equation for the diagonal elements: ```math -a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_b m_b (-d_{bi}) \nabla W_{ib}. +a_{ii} = \sum_f m_f ( d_{ii} - d_{fi}) \cdot \nabla W_{if} + \sum_b m_b d_{ii} \cdot \nabla W_{ib}. ``` From this point forward, the computation of the coefficients required for the Jacobi scheme @@ -263,11 +267,12 @@ When using pressure mirroring, the pressure value ``p_b`` of a boundary particle above is defined to be equal to the pressure of the corresponding fluid particle ``p_i``. In other words, the boundary particle "mirrors" the pressure of the fluid particle interacting with it. As a result, the coefficient that describes the influence of a particle's own -pressure value ``p_i`` ​must also include contributions from boundary particles. Therefore, -the equation for calculating the coefficient ``d_{ii}`` must be adjusted as follows: +pressure value ``p_i`` must include a doubled contribution from each boundary particle. +Therefore, ``d_{ii}`` becomes ```math -d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. +d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} + - 2\Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. ``` The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: @@ -275,30 +280,29 @@ The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: ```math \begin{align*} p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} - - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \nabla W_{if} \right. \\ -& \quad - \left. \sum_b m_b \sum_f d_{if} p_f^l \nabla W_{ib} \right). + - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \cdot \nabla W_{if} \right. \\ +& \quad - \left. \sum_b m_b \left( \sum_f d_{if} p_f^l \right) \cdot \nabla W_{ib} \right). \end{align*} ``` ### Pressure Zeroing If pressure zeroing is used instead, the pressure value of a boundary particle ``p_b`` -​is assumed to be zero. Consequently, boundary particles do not contribute to the pressure -forces acting on fluid particles. -In this case, the computation of the coefficient ``d_{ii}`` remains unchanged and is given by: +is assumed to be zero. In the linear system, this removes the boundary pressure unknowns, +but boundary particles still contribute through the ``p_i/\rho_i^2`` part of the pressure +acceleration. Therefore ``d_{ii}`` is ```math -d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if}. +d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} + - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. ``` -The equation for the relaxed Jacobi iteration remains the same as in the pressure mirroring -approach. However, the contribution from boundary particles vanishes due to their zero -pressure: +The corresponding relaxed Jacobi iteration reads ```math \begin{align*} p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} - - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \nabla W_{if} \right. \\ -& \quad - \left. \sum_b m_b \sum_j d_{if} p_f^l \nabla W_{ib} \right). + - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \cdot \nabla W_{if} \right. \\ +& \quad - \left. \sum_b m_b \left( \sum_f d_{if} p_f^l \right) \cdot \nabla W_{ib} \right). \end{align*} ``` @@ -307,8 +311,8 @@ The density calculators [`AdamiPressureExtrapolation`](@ref) and [`BernoulliPres can also be used with IISPH. When using one of these pressure extrapolation methods the calculation of the PPE is exactly the same as when using pressure zeroing. -So within the linear systems the pressure values are equal to zero (``p_b=0``) and therefore -are not considered in the calculations. Only in the pressure acceleration, the extrapolated +So within the linear system the boundary pressures are treated as zero (``p_b=0``), exactly +as in pressure zeroing. Only in the pressure acceleration, the extrapolated pressure values are used for the boundary particles. For more information on these two methods, refer to the docs for the [boundary models](@ref boundary_models). @@ -318,7 +322,7 @@ The [`PressureBoundaries`](@ref) density calculator was introduced by only be used with IISPH. In the standard IISPH method the PPE is solved only for fluid particles. The pressure values for the boundary particles are then approximated, for example by using -pressure mirroing. +pressure mirroring. With `PressureBoundaries`, however, the linear system is extended to include the boundary particles as well. This means that the pressure values of both the fluid and the boundary particles are computed directly by solving the PPE. @@ -333,22 +337,22 @@ also solved as part of the linear system. This leads to the following condition for the boundary particles ``b``: ```math -\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \nabla W_{if} + \Delta t^2 \sum_b m_b \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_b^p(t)}{m_b} \right) \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \cdot \nabla W_{if} + \Delta t^2 \sum_b m_b \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_b^p(t)}{m_b} \right) \cdot \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. ``` -Note that in this case ``i`` is a boundary particle,``f`` are its fluid neighbors, and ``b`` +Note that in this case ``i`` is a boundary particle, ``f`` are its fluid neighbors, and ``b`` its boundary neighbors. Since the pressure force for boundary particles is zero (as mentioned before), and because in this case ``i`` and ``b`` are both boundary particles, the PPE simplifies to ```math -\Delta t^2 \sum_f m_f - \frac{\bm{F}_f^p(t)}{m_f} \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. +-\Delta t^2 \sum_f m_f \frac{\bm{F}_f^p(t)}{m_f} \cdot \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. ``` If we substitute the definition of the pressure force from above, we obtain ```math -\Delta t^2 \sum_f m_f \left( \sum_k m_k \left( \frac{p_f(t)}{\rho_j^2(t)} + \frac{p_k(t)}{\rho_k^2(t)} \right) \nabla W_{fk}\right) \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \sum_k m_k \left( \frac{p_f(t)}{\rho_f^2(t)} + \frac{p_k(t)}{\rho_k^2(t)} \right) \nabla W_{fk}\right) \cdot \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. ``` where ``k`` represents all neighboring particles (fluid and boundary) of fluid particle ``f``. @@ -358,23 +362,23 @@ indirectly as neighbors of fluid particles), their ``d_{ii}`` values are zero. T diagonal elements simplify to ```math -a_{ii} = \sum_f \left( -d_{fi}\right) \nabla W_{if}. +a_{ii} = - \sum_f m_f d_{fi} \cdot \nabla W_{if}. ``` The off-diagonal term ``\sum_{j \neq i} a_{ij} p_j`` in the relaxed Jacobi iteration for boundary particles takes the following form ```math -\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( d_{ff} - \sum_{f_j} d_{ff_j}p_{f_j}\right) \nabla W_{if}. +\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( - d_{ff} p_f - \sum_{k \neq i} d_{fk}p_k \right) \cdot \nabla W_{if}. ``` But not only the addition of the boundary particles to the linear system changes when using pressure boundaries, also the PPE for the fluid particles is changing slightly. Since the boundary particles have now their own pressure values, ``p_b`` is no longer -eliminated (as in pressure zeroing or pressure extrapolation with ``p_b=0``) nor simply +eliminated (as in pressure zeroing or pressure extrapolation with ``p_b=0``) nor simply replaced by the fluid particle's pressure (as in pressure mirroring with ``p_b=p_i``). -Instead, ``p_b``remains part of the PPE. +Instead, ``p_b`` remains part of the PPE. This has no effect on the diagonal elements ``a_{ii}``, which are identical to those of the other density calculators. The ``d_{ii}`` values are also computed in the same way as for @@ -383,11 +387,11 @@ However, the off-diagonal term ``\sum_{j \neq i} a_{ij} p_j`` in the relaxed Ja changes slightly. For pressure boundaries it takes the form ```math - \sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( \sum_k d_{ik} p_k - d_{ff} p_f - \sum_{k \neq f} d_{fk} p_k \right) \nabla W_{if} \\ - + \sum_b m_b \left( \sum_k d_{ik} p_k \right) \nabla W_{ib}, +\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( \sum_k d_{ik} p_k - d_{ff} p_f - \sum_{k \neq f} d_{fk} p_k \right) \cdot \nabla W_{if} ++ \sum_b m_b \left( \sum_k d_{ik} p_k \right) \cdot \nabla W_{ib}, ``` -where ``k``represents all neighboring particles of ``i`` (both fluid and boundary). +where ``k`` represents all neighboring particles of ``i`` (both fluid and boundary). ```@docs PressureBoundaries -``` \ No newline at end of file +``` diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md index 86167d164b..56b30ebde0 100644 --- a/docs/src/systems/total_lagrangian_sph.md +++ b/docs/src/systems/total_lagrangian_sph.md @@ -106,16 +106,16 @@ Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "penalty_force ## Viscosity Another technique that is used to correct the hourglass instability is artificial viscosity. -Hereby, a viscosity term designed for fluids (see [Viscosity](@ref viscosity_sph)) is applied. -First, the force ``f_{ab}^{\text{fluid}}`` exerted by particle ``b`` on particle ``a`` +Here, a viscosity term designed for fluids (see [Viscosity](@ref viscosity_sph)) is applied. +First, the force ``\bm{F}_{ab}^{\text{fluid}}`` exerted by particle ``b`` on particle ``a`` due to artificial viscosity is computed as if both particles were fluid particles (see [Viscosity](@ref viscosity_sph) for the relevant equations). Then, according to [Lin et al. (2015)](@cite Lin2015), this force can be applied to TLSPH with the following conversion: ```math -f_{ab}^{\text{AV}} = \det(F_a) F_a^{-1} f_{ab}^{\text{fluid}}, +\bm{F}_{ab}^{\text{AV}} = \det(\bm{F}_a) \bm{F}_a^{-T} \bm{F}_{ab}^{\text{fluid}}, ``` -where ``F_a`` is the deformation gradient at particle ``a``. +where ``\bm{F}_a`` is the deformation gradient at particle ``a``. We found that artificial viscosity is not effective at correcting the incorrect particle positions due to hourglass modes. diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 2c69f4b3c9..7044ef9efa 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -52,7 +52,8 @@ pressure field. It is highly recommended to use density diffusion when using WCS ### Formulation All density diffusion terms extend the continuity equation (see [`ContinuityDensity`](@ref)) -by an additional term +by an additional term. In the literature, this is typically written for a fixed smoothing +length ``h`` as ```math \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}, @@ -62,12 +63,20 @@ the density diffusion method (see [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) for available terms). Also, ``\rho_a`` denotes the density of particle ``a`` and ``r_{ab} = r_a - r_b`` is the difference of the coordinates, ``v_{ab} = v_a - v_b`` of the velocities of particles -``a`` and ``b``. +``a`` and ``b``. When particle-wise smoothing lengths are used, the corresponding +pairwise form is +```math +\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla W_{ab} + + \delta c \sum_{b} \bar{h}_{ab} V_b \psi_{ab} \cdot \nabla W_{ab}, +``` +with ``\bar{h}_{ab} = \frac{1}{2}(h_a + h_b)``. TrixiParticles.jl uses this pairwise +average, which is the natural extension of the standard formula to variable smoothing +lengths. For fixed smoothing length, both expressions coincide. ### Numerical Results All density diffusion terms remove numerical noise in the pressure field and produce more -accurate results than weakly commpressible SPH without density diffusion. +accurate results than weakly compressible SPH without density diffusion. This can be demonstrated with dam break examples in 2D and 3D. Here, ``δ = 0.1`` has been used for all terms. Note that, due to added stability, the adaptive time integration method that was used here @@ -129,11 +138,11 @@ 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`` +After each time step, a correction term ``\delta \bm{r}_a`` is added to the position ``\bm{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} +\delta \bm{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_a W_{ab} \frac{m_b}{\rho_a + \rho_b}, ``` where: @@ -141,14 +150,27 @@ where: - ``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, +- ``\Delta x_a`` is the target particle spacing associated with particle ``a``, +- ``W(\Delta x_a)`` is the smoothing kernel evaluated at that target particle spacing, +- ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel with respect to particle ``a``, - ``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. +In TrixiParticles.jl, the same correction is applied through a shifting velocity +```math +\delta \bm{r}_a = \Delta t \, \delta \bm{v}_a, +``` +with +```math +\delta \bm{v}_a = - v_\text{max} \frac{(2h)^2}{2\Delta x} + \sum_b \left( 1 + \frac{2}{10} \left( \frac{W_{ab}}{W(\Delta x)} \right)^4 \right) + \frac{m_b}{\rho_a + \rho_b} \nabla_a W_{ab}. +``` +This corresponds to the same PST idea, but with the commonly used constants fixed to +``R = 0.2`` and ``n = 4``. The prefactor is written in a form that keeps the magnitude of +the shifting correction consistent when the smoothing-length factor changes. In particular, +``\text{CFL} \cdot \text{Ma}`` is replaced by ``\Delta t \, v_\text{max} / h``, as explained +by [Sun et al. (2018)](@cite Sun2018) on page 29, immediately above Equation 9. The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as ``\delta^+``-SPH. @@ -178,18 +200,45 @@ is a constant background pressure field. The tilde in the second term of the right-hand side indicates that the material derivative has an advection part. -The discretized form of the last term is +In the literature, the discretized form of the last term is ```math -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, ``` -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +where ``V_a`` and ``V_b`` denote the particle volumes of particles ``a`` and ``b`` respectively. Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, -the discretization is not 0th-order consistent for **non**-uniform particle distribution, -which means that there is a non-vanishing contribution only when particles are disordered. -That also means that ``p_{\text{background}}`` occurs as pre-factor to correct -the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is in the order of the reference pressure, -but it can be chosen arbitrarily large when the time-step criterion is adjusted. +the discretization is not 0th-order consistent for non-uniform particle distributions. +This means that a non-vanishing contribution appears only when the particles are disordered, +so ``p_{\text{background}}`` acts as a prefactor that regularizes the trajectories and promotes +more uniform particle distributions. + +In TrixiParticles.jl, this background-pressure contribution is evaluated through the +currently selected pressure-acceleration formulation instead of hard-coding the +specific ``(V_a^2 + V_b^2)`` discretization. This keeps the TVF term consistent with the +pressure discretization used in the rest of the scheme and automatically adapts it when a +different pressure-acceleration formulation is chosen. For the default +continuity-density formulation, +```math +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d} t}\right|_{p} += - \sum_b m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla_a W_{ab}, +``` +setting ``p_a = p_b = 1`` gives the discrete background-pressure operator +```math +- \sum_b \frac{2m_b}{\rho_a \rho_b} \nabla_a W_{ab}. +``` +In the TVF update this operator is multiplied by the background pressure and the time-step +factor from the transport-velocity correction. TrixiParticles.jl then removes the explicit +dependence on ``\Delta t`` by using the CFL estimate employed by [Adami et al. (2013)](@cite Adami2013), +```math +\Delta t \leq \frac{1}{4} \frac{h}{c_s}, +``` +as an equality, so that the resulting shifting-velocity contribution becomes +```math +\delta \bm{v}_a = - \frac{p_{\text{background}}}{8} \frac{h}{c_s} +\sum_b \frac{2m_b}{\rho_a \rho_b} \nabla_a W_{ab}, +``` +where ``h`` is the smoothing length and ``c_s`` is the speed of sound. This explains +both why the implemented pairwise factor differs from ``V_a^2 + V_b^2`` and why the +factor ``h / c_s`` appears in the final expression. The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is @@ -200,18 +249,36 @@ where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. -The discretized form of the momentum equation for a particle ``a`` reads as +The discretized form of the momentum equation for a particle ``a`` reads ```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} += \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) +\left[ -\tilde{p}_{ab} \nabla_a W_{ab} ++ \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. ``` Here, ``\tilde{p}_{ab}`` is the density-weighted pressure ```math \tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, ``` -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` -and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors -for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``, -as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. +with ``\rho_a``, ``\rho_b`` and ``p_a``, ``p_b`` denoting the densities and pressures +of particles ``a`` and ``b``, respectively. + +As for the background-pressure term above, TrixiParticles.jl evaluates the additional +convection term through the selected pressure-acceleration formulation instead of +hard-coding the ``(V_a^2 + V_b^2)`` form. For the default continuity-density +formulation, this gives +```math +\left.\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t}\right|_{\bm{A}} += - \sum_b \frac{m_b}{\rho_a \rho_b} +\left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab}. +``` +Here, ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors of particles ``a`` and ``b``, +with, for example, +```math +\bm{A}_a = \rho_a \bm{v}_a \left(\tilde{\bm{v}}_a - \bm{v}_a\right)^T. +``` +Thus, the implemented form is again the operator-consistent version of the literature +discretization for the pressure-acceleration formulation used by the scheme. To apply the TVF, use the keyword argument `shifting_technique` in the constructor of a system that supports it. @@ -247,13 +314,13 @@ Only the combination of PST and TIC is able to produce physical results. 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} +\bm{f}_{ab} = -m_a m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla_a W_{ab} ``` for the WCSPH method with [`ContinuityDensity`](@ref). -The TIC formulation changes this force to +The TIC formulation changes this term to ```math -f_{ab} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla W_{ab}. +\bm{f}_{ab}^{\mathrm{TIC}} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla_a W_{ab}. ``` Note that this formulation is asymmetric and sacrifices conservation of linear and angular momentum. diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index 006c83261d..04f68e990e 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -106,7 +106,7 @@ half step for ``v``, yielding u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, \operatorname{drift}(v^0, u^0, t^0), \\ v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\ v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\ -u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t). +u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1/2}, t^0 + \Delta t). \end{align*} ``` This scheme is implemented in OrdinaryDiffEq.jl as `LeapfrogDriftKickDrift` and yields @@ -132,7 +132,7 @@ v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\ \rho^{1/2} &= \rho^0 + \frac{1}{2} \Delta t\, R(v^0, u^0, t^0), \\ v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\ \rho^1 &= \rho^0 \frac{2 - \varepsilon^{1/2}}{2 + \varepsilon^{1/2}}, \\ -u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t), +u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1/2}, t^0 + \Delta t), \end{align*} ``` where