Mixed kernel gradient

examples/fluid/dam_break_2d.jl

@@ -16,7 +16,7 @@ fluid_density = 1000.0
 # Simulation settings
 fluid_particle_spacing = 0.02
 smoothing_length = 1.2 * fluid_particle_spacing
-boundary_layers = 3
+boundary_layers = 4
 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model
 spacing_ratio = 1
 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio

examples/fluid/dam_break_2d_corrections.jl

@@ -1,21 +1,26 @@
 using TrixiParticles
 
 fluid_density = 1000.0
+atmospheric_pressure = 100000.0
+gravity = 9.81
+
+initial_fluid_size = (2.0, 1.0)
+sound_speed = 20 * sqrt(gravity * initial_fluid_size[2])
 
 particle_spacing = 0.05
 smoothing_length = 1.15 * particle_spacing
+tspan = (0.0, 5.7 / sqrt(gravity))
 
 boundary_density_calculator = SummationDensity()
 
-gravity = 9.81
-tspan = (0.0, 5.7 / sqrt(gravity))
-
 correction_dict = Dict(
     "no_correction" => Nothing(),
     "shepard_kernel_correction" => ShepardKernelCorrection(),
     "akinci_free_surf_correction" => AkinciFreeSurfaceCorrection(fluid_density),
     "kernel_gradient_summation_correction" => KernelGradientCorrection(),
     "kernel_gradient_continuity_correction" => KernelGradientCorrection(),
+    "blended_gradient_summation_correction" => BlendedGradientCorrection(0.5),
+    "blended_gradient_continuity_correction" => BlendedGradientCorrection(0.5),
 )
 
 density_calculator_dict = Dict(
@@ -24,6 +29,18 @@ density_calculator_dict = Dict(
     "akinci_free_surf_correction" => SummationDensity(),
     "kernel_gradient_summation_correction" => SummationDensity(),
     "kernel_gradient_continuity_correction" => ContinuityDensity(),
+    "blended_gradient_summation_correction" => SummationDensity(),
+    "blended_gradient_continuity_correction" => ContinuityDensity(),
+)
+
+smoothing_kernel_dict = Dict(
+    "no_correction" => SchoenbergCubicSplineKernel{2}(),
+    "shepard_kernel_correction" => SchoenbergCubicSplineKernel{2}(),
+    "akinci_free_surf_correction" => SchoenbergCubicSplineKernel{2}(),
+    "kernel_gradient_summation_correction" => SchoenbergCubicSplineKernel{2}(),
+    "kernel_gradient_continuity_correction" => SchoenbergCubicSplineKernel{2}(),
+    "blended_gradient_summation_correction" => SchoenbergQuinticSplineKernel{2}(),
+    "blended_gradient_continuity_correction" => SchoenbergQuinticSplineKernel{2}(),
 )
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
@@ -41,6 +58,7 @@ state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pr
 for correction_name in keys(correction_dict)
     local fluid_density_calculator = density_calculator_dict[correction_name]
     local correction = correction_dict[correction_name]
+    local smoothing_kernel = smoothing_kernel_dict[correction_name]
 
     println("="^100)
     println("fluid/dam_break_2d.jl with ", correction_name)
@@ -51,6 +69,6 @@ for correction_name in keys(correction_dict)
                   boundary_density_calculator=boundary_density_calculator,
                   fluid_density_calculator=fluid_density_calculator,
                   correction=correction, use_reinit=false,
-                  state_equation=state_equation,
+                  state_equation=state_equation, smoothing_kernel=smoothing_kernel,
                   file_prefix="$(correction_name)", tspan=tspan)
 end

examples/fluid/rectangular_tank_2d.jl

@@ -80,7 +80,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # become extremely large when fluid particles are very close to boundary particles,
 # and the time integration method interprets this as an instability.
 sol = solve(ode, RDPK3SpFSAL49(),
-            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
-            reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-5, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
             save_everystep=false, callback=callbacks);

examples/fsi/dam_break_2d.jl

@@ -20,11 +20,11 @@ beta = 1
 boundary_layers = 3
 
 water_width = 0.146
-water_height = 2water_width
+water_height = 2 * water_width
 water_density = 1000.0
 
-tank_width = 4water_width
-tank_height = 4water_width
+tank_width = 4 * water_width
+tank_height = 4 * water_width
 
 sound_speed = 20 * sqrt(9.81 * water_height)
 

src/TrixiParticles.jl

@@ -52,7 +52,8 @@ export examples_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape
 export VoxelSphere, RoundSphere, reset_wall!
-export ShepardKernelCorrection, KernelGradientCorrection, AkinciFreeSurfaceCorrection
+export ShepardKernelCorrection, KernelGradientCorrection, AkinciFreeSurfaceCorrection,
+       GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles
 
 end # module

src/general/corrections.jl

@@ -1,5 +1,6 @@
-# Sorted in order of computational cost
+using LinearAlgebra
 
+# Sorted in order of computational cost
 @doc raw"""
     AkinciFreeSurfaceCorrection(rho0)
 
@@ -11,6 +12,16 @@ near free surfaces to counter this effect.
 It's important to note that this correlation is unphysical and serves as an approximation.
 The computation time added by this method is about 2-3%.
 
+Mathematically the idea is quite simple if we have a SPH particle in the middle of a volume
+at rest its density will be identical to the rest density `rho0`. If we now consider a SPH
+particle at a free surface at rest it will have neighbor missing in the direction normal to
+the surface, which will result in a lower density, if we calculate the correction factor `k`
+```math
+k = rho_0/rho_{mean}
+```
+this value will be about ~1.5 for particles at the free surface and can than be used to increase
+the pressure and viscosity accordingly.
+
 # Arguments
 - `rho0`: Reference density.
 
@@ -119,6 +130,33 @@ especially for free surfaces.
 """
 struct KernelGradientCorrection end
 
+@doc raw"""
+    MixedKernelGradientCorrection()
+
+Combines 'GradientCorrection' and 'KernelGradientCorrection' which results in a 1st order accurate SPH method.
+
+# Notes:
+- Stability issues as especially when particles separate into small clusters.
+- Doubles the computational effort.
+
+## References:
+- J. Bonet, T.-S.L. Lok.
+  "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations".
+  In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97-115.
+  [doi: 10.1016/S0045-7825(99)00051-1](https://doi.org/10.1016/S0045-7825(99)00051-1)
+- Mihai Basa, Nathan Quinlan, Martin Lastiwka.
+  "Robustness and accuracy of SPH formulations for viscous flow".
+  In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127-1148.
+  [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927)
+"""
+struct MixedKernelGradientCorrection
+    use_factorization::Bool
+
+    function MixedKernelGradientCorrection(; use_factorization=false)
+        return new{}(use_factorization)
+    end
+end
+
 function kernel_correction_coefficient(system, particle)
     return system.cache.kernel_correction_coefficient[particle]
 end
@@ -175,7 +213,8 @@ end
 
 function compute_correction_values!(system, system_index, v, u, v_ode, u_ode, semi,
                                     ::Union{SummationDensity, ContinuityDensity},
-                                    ::KernelGradientCorrection)
+                                    ::Union{KernelGradientCorrection,
+                                            MixedKernelGradientCorrection})
     (; systems, neighborhood_searches) = semi
     (; cache) = system
     (; kernel_correction_coefficient, dw_gamma) = cache
@@ -249,6 +288,8 @@ aiming to make the corrected gradient more accurate, especially near domain boun
 # Notes:
 - Stability issues as especially when particles separate into small clusters.
 - Doubles the computational effort.
+- Better stability with smoothing Kernels with larger support and Continuity e.g. 'SchoenbergQuinticSplineKernel' or 'WendlandC6Kernel'.
+- Set dt_max =< 1e-3
 
 ## References:
 - J. Bonet, T.-S.L. Lok.
@@ -260,11 +301,41 @@ aiming to make the corrected gradient more accurate, especially near domain boun
   In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127-1148.
   [doi: 10.1002/fld.1927](https://doi.org/10.1002/fld.1927)
 """
-struct GradientCorrection end
+struct GradientCorrection
+    use_factorization::Bool
+
+    function GradientCorrection(; use_factorization=false)
+        return new{}(use_factorization)
+    end
+end
+
+@doc raw"""
+    BlendedGradientCorrection()
+
+To increase stability this calculates a blended gradient, which reduces the stability issues of the 'GradientCorrection' method.
+
+This calculates the following
+```math
+\nabla^\tilde A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i
 
-function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, system,
-                                             coordinates)
-    (; mass, material_density) = system
+```
+with $\lambda$ being the blending factor.
+
+"""
+struct BlendedGradientCorrection{ELTYPE <: Real}
+    blending_factor   :: ELTYPE
+    use_factorization :: Bool
+
+    function BlendedGradientCorrection(blending_factor; use_factorization=false)
+        return new{eltype(blending_factor)}(blending_factor, use_factorization)
+    end
+end
+
+function compute_gradient_correction_matrix!(corr_matrix::AbstractArray,
+                                             neighborhood_search, system,
+                                             coordinates, density_fun;
+                                             use_factorization=false)
+    (; mass, smoothing_length, smoothing_kernel, correction) = system
 
     set_zero!(corr_matrix)
 
@@ -278,22 +349,37 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s
         # 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'
+        grad_kernel = nothing
+        if correction isa MixedKernelGradientCorrection
+            grad_kernel = corrected_kernel_grad(smoothing_kernel, pos_diff, distance,
+                                                smoothing_length,
+                                                KernelGradientCorrection(), system,
+                                                particle)
+        else
+            grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
+        end
+
+        L = volume * grad_kernel * pos_diff'
 
         @inbounds for j in 1:ndims(system), i in 1:ndims(system)
-            corr_matrix[i, j, particle] -= result[i, j]
+            corr_matrix[i, j, particle] -= L[i, j]
         end
     end
 
-    @threaded for particle in eachparticle(system)
-        L = correction_matrix(system, particle)
-        result = inv(L)
-
-        @inbounds for j in 1:ndims(system), i in 1:ndims(system)
-            corr_matrix[i, j, particle] = result[i, j]
+    if use_factorization
+        @threaded for particle in eachparticle(system)
+            L = correction_matrix(system, particle)
+            invert_factorization(corr_matrix, L, particle, system)
+        end
+    else
+        @threaded for particle in eachparticle(system)
+            L = correction_matrix(system, particle)
+            if cond(L) > 1e10
+                throw(ModelError("Correction Matrix 'L' is not well conditioned."))
+            end
+            invert(corr_matrix, L, particle, system)
         end
     end
 

src/general/density_calculators.jl

@@ -35,21 +35,6 @@ end
     return v[end, particle]
 end
 
-# *Note* that these functions are intended to internally set the density for buffer particles
-# and density correction. It cannot be used to set up an initial condition,
-# as the particle density depends on the particle positions.
-@inline function set_particle_density(particle, v, system, density)
-    set_particle_density(particle, v, system.density_calculator, system, density)
-end
-
-@inline function set_particle_density(particle, v, ::SummationDensity, system, density)
-    system.cache.density[particle] = density
-end
-
-@inline function set_particle_density(particle, v, ::ContinuityDensity, system, density)
-    v[end, particle] = density
-end
-
 function summation_density!(system, system_index, semi, u, u_ode, density;
                             particles=each_moving_particle(system))
     (; systems, neighborhood_searches) = semi

src/general/general.jl

@@ -4,6 +4,40 @@
     return du
 end
 
+@inline function invert(inverse, A, particle, system)
+    A_inv = inv(A)
+    @inbounds for j in 1:ndims(system), i in 1:ndims(system)
+        inverse[i, j, particle] = A_inv[i, j]
+    end
+end
+
+@inline function invert_factorization(inverse, A, particle, system)
+    A_factored = lu(A)
+    A_inv = similar(A)
+
+    n = ndims(system)
+    y = similar(A_inv[:, 1])
+
+    # Inverting by solving linear system for each column of the identity
+    @inbounds for j in 1:n
+        e_j = zeros(eltype(A), n)
+        e_j[j] = 1.0
+        y .= A_factored.L \ e_j
+        A_inv[:, j] .= A_factored.U \ y
+    end
+    @inbounds for j in 1:n, i in 1:n
+        inverse[i, j, particle] = A_inv[i, j]
+    end
+end
+
+struct SimulationDiverged <: Exception
+    msg::AbstractString
+end
+
+struct ModelError <: Exception
+    msg::AbstractString
+end
+
 # Note that `semidiscretization.jl` depends on the system types and has to be
 # included later.
 # `density_calculators.jl` needs to be included before `corrections.jl`.

src/general/smoothing_kernels.jl

@@ -57,6 +57,26 @@ end
            kernel_correction_coefficient(system, particle)
 end
 
+@inline function corrected_kernel_grad(kernel, pos_diff, distance, h,
+                                       corr::BlendedGradientCorrection, system, particle)
+    grad = kernel_grad(kernel, pos_diff, distance, h)
+    factor = corr.blending_factor
+    return (1 - factor) * grad + factor * correction_matrix(system, particle) * grad
+end
+
+@inline function corrected_kernel_grad(kernel, pos_diff, distance, h,
+                                       ::GradientCorrection, system, particle)
+    grad = kernel_grad(kernel, pos_diff, distance, h)
+    return correction_matrix(system, particle) * grad
+end
+
+@inline function corrected_kernel_grad(kernel, pos_diff, distance, h,
+                                       ::MixedKernelGradientCorrection, system, particle)
+    grad = corrected_kernel_grad(kernel, pos_diff, distance, h, KernelGradientCorrection(),
+                                 system, particle)
+    return correction_matrix(system, particle) * grad
+end
+
 @doc raw"""
     SchoenbergCubicSplineKernel{NDIMS}()