impl

examples/fluid/dam_break_2d_corrections.jl

@@ -1,21 +1,29 @@
 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(),
+#     "gradient_correction" => GradientCorrection(),
+# )
 
 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(),
     "gradient_correction" => GradientCorrection(),
 )
 
@@ -25,14 +33,15 @@ density_calculator_dict = Dict(
     "akinci_free_surf_correction" => SummationDensity(),
     "kernel_gradient_summation_correction" => SummationDensity(),
     "kernel_gradient_continuity_correction" => ContinuityDensity(),
+    "gradient_correction" => SummationDensity(),
 )
 
-trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-              fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length,
-              boundary_density_calculator=ContinuityDensity(),
-              fluid_density_calculator=ContinuityDensity(),
-              correction=Nothing(), use_reinit=true,
-              file_prefix="continuity_reinit", tspan=tspan)
+# trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+#               fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length,
+#               boundary_density_calculator=ContinuityDensity(),
+#               fluid_density_calculator=ContinuityDensity(),
+#               correction=Nothing(), use_reinit=true,
+#               file_prefix="continuity_reinit", tspan=tspan)
 
 # Clip negative pressure to be able to use `SummationDensity`
 state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure,

examples/fluid/rectangular_tank_2d.jl

@@ -39,7 +39,7 @@ tank = RectangularTank(particle_spacing, (water_width, water_height),
 
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
-                                             AdamiPressureExtrapolation(),
+                                             SummationDensity(),
                                              smoothing_kernel, smoothing_length)
 
 # K = 9.81 * water_height
@@ -49,11 +49,11 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
 # ==========================================================================================
 # ==== Systems
 
-density_calculator = ContinuityDensity()
+density_calculator = SummationDensity()
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, density_calculator, state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
-                                           acceleration=(0.0, gravity))
+                                           acceleration=(0.0, gravity), correction=GradientCorrection())
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
@@ -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)
-            dtmax=1e-2, # Limit stepsize to prevent crashing
+            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-4, # Limit stepsize to prevent crashing
             save_everystep=false, callback=callbacks);

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.
 
@@ -262,9 +273,9 @@ aiming to make the corrected gradient more accurate, especially near domain boun
 """
 struct GradientCorrection end
 
-function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, system,
-                                             coordinates)
-    (; mass, material_density) = system
+function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, neighborhood_search, system,
+                                             coordinates, density_fun)
+    (; mass) = system
 
     set_zero!(corr_matrix)
 
@@ -278,22 +289,48 @@ 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] / material_density[neighbor]
+        volume = mass[neighbor] / density_fun(neighbor)
+
 
         grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
-        result = volume * grad_kernel * pos_diff'
+        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)
+
+        if cond(L) > 1e10
+            #regularization
+            println("reg")
+            lambda = 1e-6
+            L = L + lambda * I
+        end
+
+
+        L_inv = inv(L)
+
+        # L_factored = lu(L)
+        # L_inv = similar(L)
+
+        # n = ndims(system)
+        # I_n = Matrix{Float64}(I, n, n) # Identity matrix of size n x n
+
+        # for j = 1:n
+        #     e_j = I_n[:, j]  # j-th column of the identity matrix
+        #     y = L_factored.L \ e_j
+        #     x_j = L_factored.U \ y
+
+        #     # Store x_j as the j-th column of L_inv
+        #     L_inv[:, j] .= x_j
+        # end
 
         @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_inv[i, j]
         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/smoothing_kernels.jl

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

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -86,12 +86,17 @@ end
 create_cache(correction, density, NDIMS, nparticles) = (;)
 
 function create_cache(::ShepardKernelCorrection, density, NDIMS, n_particles)
-    (; kernel_correction_coefficient=similar(density))
+    return (; kernel_correction_coefficient=similar(density))
 end
 
 function create_cache(::KernelGradientCorrection, density, NDIMS, n_particles)
     dw_gamma = Array{Float64}(undef, NDIMS, n_particles)
-    (; kernel_correction_coefficient=similar(density), dw_gamma)
+    return (; kernel_correction_coefficient=similar(density), dw_gamma)
+end
+
+function create_cache(::GradientCorrection, density, NDIMS, n_particles)
+    correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles)
+    return (; correction_matrix)
 end
 
 function create_cache(n_particles, ELTYPE, ::SummationDensity)
@@ -159,9 +164,13 @@ initialize!(system::WeaklyCompressibleSPHSystem, neighborhood_search) = system
 
 function update_quantities!(system::WeaklyCompressibleSPHSystem, system_index, v, u,
                             v_ode, u_ode, semi, t)
-    (; density_calculator) = system
+    (; density_calculator, correction) = system
+    (; neighborhood_searches) = semi
+
+    neighborhood_search = neighborhood_searches[system_index][system_index]
 
     compute_density!(system, system_index, u, u_ode, semi, density_calculator)
+    compute_gradient_correction_matrix!(correction, neighborhood_search, system, u, v)
 
     return system
 end
@@ -178,6 +187,20 @@ function compute_density!(system, system_index, u, u_ode, semi, ::SummationDensi
     summation_density!(system, system_index, semi, u, u_ode, density)
 end
 
+function compute_gradient_correction_matrix!(correction, neighborhood_search, system, u, v)
+    return system
+end
+
+function compute_gradient_correction_matrix!(::GradientCorrection, neighborhood_search, system, u, v)
+    (; cache) = system
+    (; correction_matrix) = cache
+
+    system_coords = current_coordinates(u, system)
+
+    compute_gradient_correction_matrix!(correction_matrix, neighborhood_search, system, system_coords, particle -> particle_density(v, system, particle))
+end
+
+
 function update_pressure!(system::WeaklyCompressibleSPHSystem, system_index, v, u,
                           v_ode, u_ode, semi, t)
     (; density_calculator, correction) = system
@@ -304,3 +327,7 @@ end
                                  system.smoothing_length,
                                  system.correction, system, particle)
 end
+
+@inline function correction_matrix(system::WeaklyCompressibleSPHSystem, particle)
+    extract_smatrix(system.cache.correction_matrix, system, particle)
+end

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -237,7 +237,7 @@ end
     return system.boundary_model.hydrodynamic_mass[particle]
 end
 
-@inline function correction_matrix(system, particle)
+@inline function correction_matrix(system::TotalLagrangianSPHSystem, particle)
     extract_smatrix(system.correction_matrix, system, particle)
 end
 @inline function deformation_gradient(system, particle)