Harden kernel gradient correction

src/general/abstract_system.jl

@@ -172,6 +172,10 @@ function update_quantities!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
+function update_density_correction!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
 function update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
@@ -180,6 +184,14 @@ function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
+function update_gradient_correction!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
+function update_surface_quantities!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
 function update_final!(system, v, u, v_ode, u_ode, semi, t; kwargs...)
     return system
 end

src/general/corrections.jl

@@ -52,7 +52,8 @@ end
     ShepardKernelCorrection()
 
 Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation
-to obtain a 0-th order accurate result, which was first proposed by [Li et al. (1996)](@cite Li1996).
+to obtain a zeroth-order consistent result (exact reproduction of constants), which was first
+proposed by [Li et al. (1996)](@cite Li1996).
 
 The kernel correction coefficient is determined by
 ```math
@@ -61,7 +62,10 @@ c(x) = \sum_{b=1} V_b W_b(x),
 where ``V_b = m_b / \rho_b`` is the volume of particle ``b``.
 
 This correction is applied with [`SummationDensity`](@ref) to correct the density and leads
-to an improvement, especially at free surfaces.
+to an improvement, especially at free surfaces. With summation density, the current one-pass
+implementation uses the provisional density in ``V_b`` and therefore reduces the free-surface
+error without guaranteeing convergence. [`DensityReinitializationCallback`](@ref) instead uses
+the independently evolved continuity density and realizes the consistent Shepard operator.
 
 !!! note
     - It is also referred to as "0th order correction".
@@ -73,7 +77,8 @@ struct ShepardKernelCorrection end
     KernelCorrection()
 
 Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation
-to obtain a 0-th order accurate result, which was first proposed by Li et al.
+to obtain a zeroth-order consistent kernel gradient (an exact zero gradient for constants),
+which was first proposed by Li et al.
 This can be further extended to obtain a kernel corrected gradient as shown by [Basa et al. (2008)](@cite Basa2008).
 
 The kernel correction coefficient is determined by
@@ -108,6 +113,16 @@ which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1
 """
 struct MixedKernelGradientCorrection end
 
+correction_density(::Any) = nothing
+correction_density(correction::ShepardKernelCorrection) = correction
+
+correction_gradient(::Nothing) = nothing
+correction_gradient(::ShepardKernelCorrection) = nothing
+correction_gradient(::AkinciFreeSurfaceCorrection) = nothing
+correction_gradient(correction) = correction
+
+correction_force(correction) = correction
+
 function kernel_correction_coefficient(system::AbstractFluidSystem, particle)
     return system.cache.kernel_correction_coefficient[particle]
 end
@@ -166,9 +181,24 @@ function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi,
         end
     end
 
+    sanitize_kernel_correction_coefficient!(kernel_correction_coefficient, system, semi)
+
     return kernel_correction_coefficient
 end
 
+function sanitize_kernel_correction_coefficient!(coefficient, system, semi)
+    minimum_coefficient = sqrt(eps(eltype(coefficient)))
+
+    @threaded semi for particle in eachparticle(system)
+        value = coefficient[particle]
+        if !isfinite(value) || value <= minimum_coefficient
+            coefficient[particle] = one(value)
+        end
+    end
+
+    return coefficient
+end
+
 function dw_gamma(system::AbstractFluidSystem, particle)
     return extract_svector(system.cache.dw_gamma, system, particle)
 end
@@ -255,9 +285,22 @@ function compute_correction_values!(system,
         end
     end
 
-    for particle in eachparticle(system), i in axes(dw_gamma, 1)
-        dw_gamma[i, particle] /= kernel_correction_coefficient[particle]
+    minimum_coefficient = sqrt(eps(eltype(kernel_correction_coefficient)))
+    @threaded semi for particle in eachparticle(system)
+        coefficient = kernel_correction_coefficient[particle]
+        if !isfinite(coefficient) || coefficient <= minimum_coefficient
+            kernel_correction_coefficient[particle] = one(coefficient)
+            for i in axes(dw_gamma, 1)
+                dw_gamma[i, particle] = zero(eltype(dw_gamma))
+            end
+        else
+            for i in axes(dw_gamma, 1)
+                dw_gamma[i, particle] /= coefficient
+            end
+        end
     end
+
+    return kernel_correction_coefficient
 end
 
 @doc raw"""

src/general/semidiscretization.jl

@@ -633,7 +633,13 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
 
     update_implicit_sph!(semi, v_ode, u_ode, t)
 
-    # Perform correction and pressure calculation
+    # Correction moments can use densities from every interacting system, so density
+    # correction has to be a global phase.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_density_correction!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
+    # Fluid pressure must be available before boundary pressure interpolation.
     foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
         update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     end
@@ -644,6 +650,18 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
         update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     end
 
+    # Boundary interpolation can update boundary density. Assemble all gradient corrections
+    # only after every interacting system exposes its final density.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_gradient_correction!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
+    # Surface quantities can depend on corrected gradients and must be complete for every
+    # system before curvature and stress are computed in `update_final!`.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_surface_quantities!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
     # Final update step for all remaining systems
     foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
         update_final!(system, v, u, v_ode, u_ode, semi, t)

src/schemes/boundary/wall_boundary/dummy_particles.jl

@@ -237,7 +237,7 @@ function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particl
 end
 
 function create_cache_model(::KernelCorrection, density, NDIMS, n_particles)
-    dw_gamma = Array{Float64}(undef, NDIMS, n_particles)
+    dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles)
     return (; kernel_correction_coefficient=similar(density), dw_gamma)
 end
 
@@ -393,21 +393,62 @@ end
 
 @inline function update_pressure!(boundary_model::BoundaryModelDummyParticles,
                                   system, v, u, v_ode, u_ode, semi)
-    (; correction, density_calculator) = boundary_model
+    (; density_calculator) = boundary_model
 
     compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
 
-    # These are only computed when using corrections
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
-    compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode,
-                                        semi)
-    # `kernel_correct_density!` only performed for `SummationDensity`
-    kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction,
+    return boundary_model
+end
+
+@inline function update_density_correction!(boundary_model::BoundaryModelDummyParticles,
+                                            system, v, u, v_ode, u_ode, semi)
+    (; correction, density_calculator) = boundary_model
+    density_correction = correction_density(correction)
+
+    compute_boundary_correction_values!(boundary_model, system, density_correction, u,
+                                        v_ode, u_ode, semi)
+    kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi,
+                            density_correction,
                             density_calculator)
 
     return boundary_model
 end
 
+@inline function update_gradient_correction!(boundary_model::BoundaryModelDummyParticles,
+                                             system, v, u, v_ode, u_ode, semi)
+    gradient_correction = correction_gradient(boundary_model.correction)
+
+    compute_boundary_correction_values!(boundary_model, system, gradient_correction, u,
+                                        v_ode, u_ode, semi)
+    compute_gradient_correction_matrix!(gradient_correction, boundary_model, system, u,
+                                        v_ode, u_ode, semi)
+
+    return boundary_model
+end
+
+@inline function compute_boundary_correction_values!(boundary_model, system, correction, u,
+                                                     v_ode, u_ode, semi)
+    return boundary_model
+end
+
+function compute_boundary_correction_values!(boundary_model, system,
+                                             ::ShepardKernelCorrection, u,
+                                             v_ode, u_ode, semi)
+    return compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode,
+                                  semi,
+                                  boundary_model.cache.kernel_correction_coefficient)
+end
+
+function compute_boundary_correction_values!(boundary_model, system,
+                                             correction::Union{KernelCorrection,
+                                                               MixedKernelGradientCorrection},
+                                             u, v_ode, u_ode, semi)
+    return compute_correction_values!(system, correction, current_coordinates(u, system),
+                                      v_ode, u_ode, semi,
+                                      boundary_model.cache.kernel_correction_coefficient,
+                                      boundary_model.cache.dw_gamma)
+end
+
 function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi,
                                  correction, density_calculator)
     return boundary_model
@@ -428,13 +469,13 @@ function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
                                                          MixedKernelGradientCorrection},
                                              boundary_model,
                                              system, u, v_ode, u_ode, semi)
-    (; cache, correction, smoothing_kernel) = boundary_model
+    (; cache, smoothing_kernel) = boundary_model
     (; correction_matrix) = cache
 
     system_coords = current_coordinates(u, system)
 
     compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
-                                        v_ode, u_ode, semi, correction, smoothing_kernel)
+                                        v_ode, u_ode, semi, corr, smoothing_kernel)
 end
 
 function compute_density!(boundary_model, ::SummationDensity, system, v, u, v_ode, u_ode,

src/schemes/boundary/wall_boundary/system.jl

@@ -218,6 +218,13 @@ function update_quantities!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi
     return system
 end
 
+function update_density_correction!(system::WallBoundarySystem{<:BoundaryModelDummyParticles},
+                                    v, u, v_ode, u_ode, semi, t)
+    update_density_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 # This update depends on the computed quantities of the fluid system and therefore
 # has to be in `update_boundary_interpolation!` after `update_quantities!`.
 function update_boundary_interpolation!(system::WallBoundarySystem, v, u, v_ode, u_ode,
@@ -231,6 +238,13 @@ function update_boundary_interpolation!(system::WallBoundarySystem, v, u, v_ode,
     return system
 end
 
+function update_gradient_correction!(system::WallBoundarySystem{<:BoundaryModelDummyParticles},
+                                     v, u, v_ode, u_ode, semi, t)
+    update_gradient_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 function write_u0!(u0, ::WallBoundarySystem)
     return u0
 end
@@ -329,7 +343,7 @@ function system_smoothing_kernel(system::WallBoundarySystem{<:BoundaryModelDummy
 end
 
 function system_correction(system::WallBoundarySystem{<:BoundaryModelDummyParticles})
-    return system.boundary_model.correction
+    return correction_gradient(system.boundary_model.correction)
 end
 
 @inline function density_calculator(system::WallBoundarySystem)

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -4,6 +4,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                    particle_system::EntropicallyDampedSPHSystem,
                    neighbor_system, semi)
     (; sound_speed, density_calculator, correction, nu_edac) = particle_system
+    gradient_correction = correction_gradient(correction)
 
     system_coords = current_coordinates(u_particle_system, particle_system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
@@ -63,7 +64,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             particle, neighbor,
                                             m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a,
                                             rho_b, pos_diff, distance, grad_kernel,
-                                            correction)
+                                            gradient_correction)
 
         dv_particle = Ref(dv_pressure)
         @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
@@ -77,7 +78,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                particle_system, neighbor_system,
                                v_particle_system, v_neighbor_system,
                                particle, neighbor, m_a, m_b, rho_a, rho_b, v_a, v_b,
-                               pos_diff, distance, grad_kernel, correction)
+                               pos_diff, distance, grad_kernel, gradient_correction)
 
         @inbounds surface_tension_force!(dv_particle, surface_tension_a,
                                          surface_tension_b,