Mixed kernel gradient

examples/fluid/dam_break_2d.jl

@@ -46,7 +46,7 @@ smoothing_kernel = WendlandC2Kernel{2}()
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
-#density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
+# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,

examples/fluid/dam_break_2d_corrections.jl

@@ -7,7 +7,7 @@ particle_spacing = 0.05
 tspan = (0.0, 5.7 / sqrt(9.81))
 
 correction_dict = Dict(
-    "no_correction" => Nothing(),
+    "no_correction" => nothing,
     "shepard_kernel_correction" => ShepardKernelCorrection(),
     "akinci_free_surf_correction" => AkinciFreeSurfaceCorrection(fluid_density),
     "kernel_gradient_summation_correction" => KernelCorrection(),

src/general/corrections.jl

@@ -361,7 +361,7 @@ struct BlendedGradientCorrection{ELTYPE <: Real}
     end
 end
 
-# called by DensityDiffusion and TLSPH
+# Called only by DensityDiffusion and TLSPH
 function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search,
                                              system, coordinates, density_fun)
     (; mass) = system
@@ -393,7 +393,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search,
 end
 
 function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
-                                             coordinates, u_ode, v_ode, semi,
+                                             coordinates, v_ode, u_ode, semi,
                                              correction, smoothing_length, smoothing_kernel)
     set_zero!(corr_matrix)
 
@@ -448,20 +448,6 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
 end
 
 function correction_matrix_inversion_step!(corr_matrix, system)
-    @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 pseudo_invert(inverse, A, particle, system)
-        A_inv = pinv(A)
-        @inbounds for j in 1:ndims(system), i in 1:ndims(system)
-            inverse[i, j, particle] = A_inv[i, j]
-        end
-    end
-
     @threaded for particle in eachparticle(system)
         L = extract_smatrix(corr_matrix, system, particle)
         norm_ = norm(L)
@@ -481,9 +467,15 @@ function correction_matrix_inversion_step!(corr_matrix, system)
 
         det_ = abs(det(L))
         @fastmath if det_ < 1e-6 * norm_
-            pseudo_invert(corr_matrix, L, particle, system)
+            L_inv = pinv(L)
+            @inbounds for j in 1:ndims(system), i in 1:ndims(system)
+                corr_matrix[i, j, particle] = L_inv[i, j]
+            end
             continue
         end
-        invert(corr_matrix, L, particle, system)
+        L_inv = inv(L)
+        @inbounds for j in 1:ndims(system), i in 1:ndims(system)
+            corr_matrix[i, j, particle] = L_inv[i, j]
+        end
     end
 end

src/general/general.jl

@@ -1,12 +1,9 @@
 abstract type System{NDIMS} end
 
-# WCSPH EDAC
 abstract type FluidSystem{NDIMS} <: System{NDIMS} end
 
-# TLSPH
 abstract type SolidSystem{NDIMS} <: System{NDIMS} end
 
-# all Boundary Condition Systems
 abstract type BoundarySystem{NDIMS} <: System{NDIMS} end
 
 timer_name(::FluidSystem) = "fluid"

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -308,12 +308,12 @@ end
                                                              GradientCorrection,
                                                              BlendedGradientCorrection,
                                                              MixedKernelGradientCorrection})
+    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
 
     # Use `p_a` as pressure for both particles with `PressureMirroring`
     return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_a,
                                             rho_a, rho_b, pos_diff, distance, grad_kernel,
-                                            particle_system, neighbor, neighbor_system,
-                                            fluid_density_calculator)
+                                            W_b, fluid_density_calculator)
 end
 
 @inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
@@ -326,10 +326,11 @@ end
                                                              GradientCorrection,
                                                              BlendedGradientCorrection,
                                                              MixedKernelGradientCorrection})
+    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
+
     return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b,
                                             rho_a, rho_b, pos_diff, distance, grad_kernel,
-                                            particle_system, neighbor, neighbor_system,
-                                            fluid_density_calculator)
+                                            W_b, fluid_density_calculator)
 end
 
 @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle)
@@ -375,8 +376,8 @@ end
     compute_correction_values!(system,
                                correction, v, u, v_ode, u_ode, semi, density_calculator)
 
-    compute_gradient_correction_matrix!(correction, boundary_model, system, u, u_ode,
-                                        v_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,
@@ -398,22 +399,22 @@ function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi,
 end
 
 function compute_gradient_correction_matrix!(correction, boundary_model, system, u,
-                                             u_ode, v_ode, semi)
+                                             v_ode, u_ode, semi)
     return system
 end
 
 function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
                                                          BlendedGradientCorrection,
                                                          MixedKernelGradientCorrection},
                                              boundary_model,
-                                             system, u, u_ode, v_ode, semi)
+                                             system, u, v_ode, u_ode, semi)
     (; cache, correction, smoothing_kernel, smoothing_length) = boundary_model
     (; correction_matrix) = cache
 
     system_coords = current_coordinates(u, system)
 
     compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
-                                        u_ode, v_ode, semi, correction, smoothing_length,
+                                        v_ode, u_ode, semi, correction, smoothing_length,
                                         smoothing_kernel)
 end
 

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -77,7 +77,8 @@ end
                                        density_calculator,
                                        correction)
 
-    # By default, just call the pressure acceleration formulation corresponding to the density calculator
+    # Without correction, the kernel gradient is symmetric, so call the symmetric
+    # pressure acceleration formulation corresponding to the density calculator.
     return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, rho_b,
                                            grad_kernel, density_calculator)
 end
@@ -91,22 +92,24 @@ end
                                                          GradientCorrection,
                                                          BlendedGradientCorrection,
                                                          MixedKernelGradientCorrection})
+    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
 
-    # By default, just call the pressure acceleration formulation corresponding to a locally corrected gradient
+    # With correction, the kernel gradient is not necessarily symmetric, so call the
+    # asymmetric pressure acceleration formulation corresponding to the density calculator.
     return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b,
                                             rho_a, rho_b, pos_diff, distance, grad_kernel,
-                                            particle_system, neighbor,
-                                            neighbor_system, density_calculator)
+                                            W_b, density_calculator)
 end
 
 # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
 # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
 # used with `ContinuityDensity` with the formulation `\rho_a * \sum m_b / \rho_b ...`.
 # This can also be seen in the tests for total energy conservation, which fail with the
 # other `pressure_acceleration` form.
+# We assume symmetry of the kernel gradient in this formulation. See below for the
+# asymmetric version.
 @inline function pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a,
-                                                 rho_b,
-                                                 grad_kernel, ::ContinuityDensity)
+                                                 rho_b, grad_kernel, ::ContinuityDensity)
     return (-m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel) * pressure_correction
 end
 
@@ -115,28 +118,27 @@ end
 # used with `SummationDensity`.
 # This can also be seen in the tests for total energy conservation, which fail with the
 # other `pressure_acceleration` form.
+# We assume symmetry of the kernel gradient in this formulation. See below for the
+# asymmetric version.
 @inline function pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a,
-                                                 rho_b,
-                                                 grad_kernel, ::SummationDensity)
+                                                 rho_b, grad_kernel, ::SummationDensity)
     return (-m_b * (p_a / rho_a^2 + p_b / rho_b^2) * grad_kernel) * pressure_correction
 end
 
+# Same as above, but not assuming symmetry of the kernel gradient. To be used with
+# corrections that do not produce a symmetric kernel gradient.
 @inline function pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b,
                                                   rho_a, rho_b, pos_diff, distance, W_a,
-                                                  particle_system, neighbor,
-                                                  neighbor_system, ::SummationDensity)
-    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
-
-    return (-m_b * (p_a / rho_a^2 * W_a - p_b / rho_b^2 * W_b)) * pressure_correction
+                                                  W_b, ::ContinuityDensity)
+    return -m_b / (rho_a * rho_b) * (p_a * W_a - p_b * W_b) * pressure_correction
 end
 
+# Same as above, but not assuming symmetry of the kernel gradient. To be used with
+# corrections that do not produce a symmetric kernel gradient.
 @inline function pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b,
                                                   rho_a, rho_b, pos_diff, distance, W_a,
-                                                  particle_system, neighbor,
-                                                  neighbor_system, ::ContinuityDensity)
-    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
-
-    return -m_b / (rho_a * rho_b) * (p_a * W_a - p_b * W_b) * pressure_correction
+                                                  W_b, ::SummationDensity)
+    return (-m_b * (p_a / rho_a^2 * W_a - p_b / rho_b^2 * W_b)) * pressure_correction
 end
 
 # With 'SummationDensity', density is calculated in wcsph/system.jl:compute_density!

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -203,7 +203,7 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od
     compute_correction_values!(system,
                                correction, v, u, v_ode, u_ode, semi, density_calculator)
 
-    compute_gradient_correction_matrix!(correction, system, u, u_ode, v_ode, semi)
+    compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi)
 
     # `kernel_correct_density!` only performed for `SummationDensity`
     kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction,
@@ -225,22 +225,22 @@ end
 
 function compute_gradient_correction_matrix!(correction,
                                              system::WeaklyCompressibleSPHSystem, u,
-                                             u_ode, v_ode, semi)
+                                             v_ode, u_ode, semi)
     return system
 end
 
 function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
                                                          BlendedGradientCorrection,
                                                          MixedKernelGradientCorrection},
                                              system::WeaklyCompressibleSPHSystem, u,
-                                             u_ode, v_ode, semi)
+                                             v_ode, u_ode, semi)
     (; cache, correction, smoothing_kernel, smoothing_length) = system
     (; correction_matrix) = cache
 
     system_coords = current_coordinates(u, system)
 
     compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
-                                        u_ode, v_ode, semi, correction, smoothing_length,
+                                        v_ode, u_ode, semi, correction, smoothing_length,
                                         smoothing_kernel)
 end
 

src/schemes/solid/total_lagrangian_sph/rhs.jl

@@ -91,10 +91,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance)
 
         # use `m_a` to get the same viscosity as for the fluid-solid direction.
-        dv_viscosity = viscosity(neighbor_system, particle_system,
-                                 v_neighbor_system, v_particle_system,
-                                 neighbor, particle,
-                                 pos_diff, distance, sound_speed, m_b, m_a, rho_mean)
+        dv_viscosity = viscosity(neighbor_system, particle_system, v_neighbor_system,
+                                 v_particle_system, neighbor, particle, pos_diff, distance,
+                                 sound_speed, m_b, m_a, rho_mean)
 
         # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles
         # corresponding to the chosen boundary model.
@@ -104,9 +103,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         # Boundary forces
         # Note: neighbor and particle pressure are switched in this call
         #       and `pressure_correction` is set to `1.0` (no correction)
-        dv_boundary = pressure_acceleration(1.0, m_b, p_b, p_a,
+        dv_boundary = pressure_acceleration(1.0, m_a, p_b, p_a,
                                             rho_b, rho_a, pos_diff, distance, grad_kernel,
-                                            particle_system, neighbor, neighbor_system,
+                                            neighbor_system, neighbor, particle_system,
                                             density_calculator, correction)
         dv_particle = dv_boundary + dv_viscosity
 

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -129,22 +129,14 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem
                       ((1 + poisson_ratio) * (1 - 2 * poisson_ratio))
         lame_mu = 0.5 * young_modulus / (1 + poisson_ratio)
 
-        return new{typeof(boundary_model),
-                   NDIMS, ELTYPE,
+        return new{typeof(boundary_model), NDIMS, ELTYPE,
                    typeof(smoothing_kernel),
-                   typeof(penalty_force)}(initial_condition,
-                                          initial_coordinates,
-                                          current_coordinates, mass,
-                                          correction_matrix,
-                                          pk1_corrected,
-                                          deformation_grad,
-                                          material_density,
-                                          n_moving_particles,
-                                          young_modulus, poisson_ratio,
-                                          lame_lambda, lame_mu,
-                                          smoothing_kernel,
-                                          smoothing_length,
-                                          acceleration_, boundary_model,
+                   typeof(penalty_force)}(initial_condition, initial_coordinates,
+                                          current_coordinates, mass, correction_matrix,
+                                          pk1_corrected, deformation_grad, material_density,
+                                          n_moving_particles, young_modulus, poisson_ratio,
+                                          lame_lambda, lame_mu, smoothing_kernel,
+                                          smoothing_length, acceleration_, boundary_model,
                                           penalty_force)
     end
 end