Implement density diffusion (#272)

* Implement density diffusion * Improve `show` * Add comment referencing a new issue * Fix `show` * Add docs * [skip ci] Consistency * Fix tests * [skip ci] Fix example tests * Fix typos in docs * Revise docs * Fix merge * Make `delta` a kwarg * Update docs * Update estimated overhead by Antuono DD * Make Molteni & Colagrossi DDT default for dam breaks * Fix tests * Implement suggestions * Reformat
trixi-framework · Nov 16, 2023 · b0b7742 · b0b7742
1 parent 5fe4af0
commit b0b7742
Show file tree

Hide file tree

Showing 17 changed files with 402 additions and 32 deletions.
diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md
@@ -25,6 +25,78 @@ Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "viscosity.jl")]
 ```
 
+## Density Diffusion
+
+Density diffusion can be used with [`ContinuityDensity`](@ref) to remove the noise in the
+pressure field. It is highly recommended to use density diffusion when using WCSPH.
+
+### Formulation
+
+All density diffusion terms extend the continuity equation (see [`ContinuityDensity`](@ref))
+by an additional term
+```math
+\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla_{r_a} W(\Vert r_{ab} \Vert, h)
+    + \delta h c \sum_{b} V_b \psi_{ab} \cdot \nabla_{r_a} W(\Vert r_{ab} \Vert, h),
+```
+where ``V_b = m_b / \rho_b`` is the volume of particle ``b`` and ``\psi_{ab}`` depends on
+the density diffusion method (see [`DensityDiffusion`](@ref) for available terms).
+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``.
+
+### 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.
+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
+can choose higher time steps in the simulations with density diffusion.
+For the cheap [`DensityDiffusionMolteniColagrossi`](@ref), this results in reduced runtime.
+
+```@raw html
+<figure>
+  <img src="https://lh3.googleusercontent.com/drive-viewer/AK7aPaBL-tqW6p9ry3NHvNnHVNufRfh_NSz0Le4vJ4n2rS-10Vr3Dkm2Cjb4T861vk6yhnvqMgS_PLXeZsNoVepIfYgpw-hlgQ=s1600" alt="density_diffusion_2d"/>
+  <figcaption>Dam break in 2D with different density diffusion terms</figcaption>
+</figure>
+```
+
+```@raw html
+<figure>
+  <img src="https://lh3.googleusercontent.com/drive-viewer/AK7aPaDKc0DCJfFH606zWFkjutMYzs70Y4Ot_33avjcIRxV3xNbrX1gqx6EpeSmysai338aRsOoqJ8B1idUs5U30SA_o12OQ=s1600" alt="density_diffusion_3d"/>
+  <figcaption>Dam break in 3D with different density diffusion terms</figcaption>
+</figure>
+```
+
+The simpler terms [`DensityDiffusionMolteniColagrossi`](@ref) and
+[`DensityDiffusionFerrari`](@ref) do not solve the hydrostatic problem and lead to incorrect
+solutions in long-running steady-state hydrostatic simulations with free surfaces
+(Antuono et al., 2012). This can be seen when running the simple rectangular tank example
+until ``t = 40`` (again using ``δ = 0.1``):
+
+```@raw html
+<figure>
+  <img src="https://lh3.googleusercontent.com/drive-viewer/AK7aPaCf1gDlbxkQjxpyffPJ-ijx-DdVxlwUVb_DLYIW4X5E0hkDeJcuAqCae6y4eDydgTKe752zWa08tKVL5yhB-ad8Uh8J=s1600" alt="density_diffusion_tank"/>
+  <figcaption>Tank in rest under gravity in 3D with different density diffusion terms</figcaption>
+</figure>
+```
+
+[`DensityDiffusionAntuono`](@ref) adds a correction term to solve this problem, but this
+term is very expensive and adds about 40--50% of computational cost.
+
+### References
+- M. Antuono, A. Colagrossi, S. Marrone.
+  "Numerical Diffusive Terms in Weakly-Compressible SPH Schemes."
+  In: Computer Physics Communications 183.12 (2012), pages 2570--2580.
+  [doi: 10.1016/j.cpc.2012.07.006](https://doi.org/10.1016/j.cpc.2012.07.006)
+
+### API
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "density_diffusion.jl")]
+```
+
 ## Corrections
 
 ```@autodocs

diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -44,10 +44,12 @@ smoothing_kernel = WendlandC2Kernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
+density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
                                            smoothing_length, viscosity=viscosity,
+                                           density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity), correction=nothing)
 
 # ==========================================================================================

diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl
@@ -6,7 +6,7 @@ using OrdinaryDiffEq
 fluid_particle_spacing = 0.08
 
 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model
-boundary_layers = 3
+boundary_layers = 4
 spacing_ratio = 1
 
 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio
@@ -32,15 +32,17 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl
 
 # ==========================================================================================
 # ==== Fluid
-smoothing_length = 1.2 * fluid_particle_spacing
-smoothing_kernel = SchoenbergCubicSplineKernel{3}()
+smoothing_length = 3.0 * fluid_particle_spacing
+smoothing_kernel = WendlandC2Kernel{3}()
 
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
+density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
                                            smoothing_length, viscosity=viscosity,
+                                           density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity, 0.0))
 
 # ==========================================================================================

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -45,6 +45,8 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel
 export StateEquationIdealGas, StateEquationCole
 export ArtificialViscosityMonaghan, ViscosityAdami
+export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
+       DensityDiffusionAntuono
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
        PressureMirroring, PressureZeroing
 export BoundaryMovement

diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl
@@ -204,6 +204,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "save final solution" => solution_saving.save_final_solution ? "yes" :
                                      "no",
             "output directory" => abspath(normpath(solution_saving.output_directory)),
+            "prefix" => solution_saving.prefix,
         ]
         summary_box(io, "SolutionSavingCallback", setup)
     end
@@ -228,6 +229,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "save final solution" => solution_saving.save_final_solution ? "yes" :
                                      "no",
             "output directory" => abspath(normpath(solution_saving.output_directory)),
+            "prefix" => solution_saving.prefix,
         ]
         summary_box(io, "SolutionSavingCallback", setup)
     end

diff --git a/src/general/corrections.jl b/src/general/corrections.jl
@@ -264,8 +264,8 @@ the correction matrix $\bm{L}_a$ is evaluated explicitly as
 struct GradientCorrection end
 
 function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, system,
-                                             coordinates)
-    (; mass, material_density) = system
+                                             coordinates, density_fun)
+    (; mass) = system
 
     set_zero!(corr_matrix)
 
@@ -274,12 +274,11 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s
                           coordinates, coordinates,
                           neighborhood_search;
                           particles=eachparticle(system)) do particle, neighbor,
-                                                             pos_diff,
-                                                             distance
+                                                             pos_diff, distance
         # Only consider particles with a distance > 0.
         distance < sqrt(eps()) && return
 
-        volume = mass[neighbor] / material_density[neighbor]
+        volume = mass[neighbor] / density_fun(neighbor)
 
         grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
         result = volume * grad_kernel * pos_diff'
@@ -290,9 +289,15 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s
     end
 
     @threaded for particle in eachparticle(system)
-        L = correction_matrix(system, particle)
+        L = extract_smatrix(corr_matrix, system, particle)
         result = inv(L)
 
+        if any(isinf.(result)) || any(isnan.(result))
+            # TODO How do we handle singularities correctly?
+            # See https://github.com/trixi-framework/TrixiParticles.jl/issues/273
+            result = one(L)
+        end
+
         @inbounds for j in 1:ndims(system), i in 1:ndims(system)
             corr_matrix[i, j, particle] = result[i, j]
         end

diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl
@@ -17,9 +17,9 @@ Density calculator to integrate the density from the continuity equation
 ```math
 \frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla_{r_a} W(\Vert r_a - r_b \Vert, h),
 ```
-where ``\rho_a`` denotes the density of particle ``a``, ``r_a`` and ``r_b`` denote the coordinates
-of particles ``a`` and ``b`` respectively, and ``v_{ab} = v_a - v_b`` is the difference of the
-velocities of particles ``a`` and ``b``.
+where ``\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``.
 """
 struct ContinuityDensity end
 

diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl
@@ -4,7 +4,7 @@
 An abstract supertype of all smoothing kernels. The type parameter `NDIMS` encodes the number
 of dimensions.
 
-We provide the following smoothing kernels:
+Currently the following smoothing kernels are available:
 
 | Smoothing Kernel                          | Compact Support   | Typ. Smoothing Length | Recommended Application | Stability |
 | :---------------------------------------- | :---------------- | :-------------------- | :---------------------- | :-------- |