diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 6a2551246..4f3b7d5fa 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -13,7 +13,7 @@ using OrdinaryDiffEq fluid_particle_spacing = 0.02 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model -boundary_layers = 4 +boundary_layers = 3 spacing_ratio = 1 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio @@ -31,7 +31,8 @@ fluid_density = 1000.0 atmospheric_pressure = 100000.0 sound_speed = 20 * sqrt(gravity * initial_fluid_size[2]) state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure, - background_pressure=atmospheric_pressure) + background_pressure=atmospheric_pressure, + clip_negative_pressure=false) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, @@ -45,6 +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) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, @@ -58,7 +60,8 @@ boundary_density_calculator = AdamiPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length, + correction=nothing) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) @@ -68,7 +71,7 @@ semi = Semidiscretization(fluid_system, boundary_system, neighborhood_search=GridNeighborhoodSearch) ode = semidiscretize(semi, tspan) -info_callback = InfoCallback(interval=100) +info_callback = InfoCallback(interval=250) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") use_reinit = false @@ -87,5 +90,5 @@ callbacks = CallbackSet(info_callback, saving_callback, density_reinit_cb) sol = solve(ode, RDPK3SpFSAL49(), 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 + dtmax=1e-3, # Limit stepsize to prevent crashing save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_2d_corrections.jl b/examples/fluid/dam_break_2d_corrections.jl index ad697d454..bd2b790af 100644 --- a/examples/fluid/dam_break_2d_corrections.jl +++ b/examples/fluid/dam_break_2d_corrections.jl @@ -3,19 +3,35 @@ using TrixiParticles fluid_density = 1000.0 particle_spacing = 0.05 -smoothing_length = 1.15 * particle_spacing -boundary_density_calculator = SummationDensity() - -gravity = 9.81 -tspan = (0.0, 5.7 / sqrt(gravity)) +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" => KernelGradientCorrection(), - "kernel_gradient_continuity_correction" => KernelGradientCorrection(), + "kernel_gradient_summation_correction" => KernelCorrection(), + "kernel_gradient_continuity_correction" => KernelCorrection(), + "blended_gradient_summation_correction" => BlendedGradientCorrection(0.5), + "blended_gradient_continuity_correction" => BlendedGradientCorrection(0.2), + "gradient_summation_correction" => GradientCorrection(), + "mixed_kernel_gradient_summation_correction" => MixedKernelGradientCorrection(), + # "gradient_continuity_correction" => GradientCorrection(), + # "mixed_kernel_gradient_continuity_correction" => MixedKernelGradientCorrection(), +) + +smoothing_length_dict = Dict( + "no_correction" => 3.0 * particle_spacing, + "shepard_kernel_correction" => 3.0 * particle_spacing, + "akinci_free_surf_correction" => 3.0 * particle_spacing, + "kernel_gradient_summation_correction" => 4.0 * particle_spacing, + "kernel_gradient_continuity_correction" => 3.5 * particle_spacing, + "blended_gradient_summation_correction" => 3.0 * particle_spacing, + "blended_gradient_continuity_correction" => 4.0 * particle_spacing, + "gradient_summation_correction" => 3.5 * particle_spacing, + "mixed_kernel_gradient_summation_correction" => 3.5 * particle_spacing, + "gradient_continuity_correction" => 4.5 * particle_spacing, + "mixed_kernel_gradient_continuity_correction" => 4.0 * particle_spacing, ) density_calculator_dict = Dict( @@ -24,23 +40,42 @@ 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(), + "gradient_summation_correction" => SummationDensity(), + "gradient_continuity_correction" => ContinuityDensity(), + "mixed_kernel_gradient_summation_correction" => SummationDensity(), + "mixed_kernel_gradient_continuity_correction" => ContinuityDensity(), +) + +smoothing_kernel_dict = Dict( + "no_correction" => WendlandC2Kernel{2}(), + "shepard_kernel_correction" => WendlandC2Kernel{2}(), + "akinci_free_surf_correction" => WendlandC2Kernel{2}(), + "kernel_gradient_summation_correction" => WendlandC6Kernel{2}(), + "kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(), + "blended_gradient_summation_correction" => WendlandC2Kernel{2}(), + "blended_gradient_continuity_correction" => WendlandC6Kernel{2}(), + "gradient_summation_correction" => WendlandC6Kernel{2}(), + "gradient_continuity_correction" => WendlandC6Kernel{2}(), + "mixed_kernel_gradient_summation_correction" => WendlandC6Kernel{2}(), + "mixed_kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(), ) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), - fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length, + fluid_particle_spacing=particle_spacing, + smoothing_length=3.0 * particle_spacing, boundary_density_calculator=ContinuityDensity(), fluid_density_calculator=ContinuityDensity(), - correction=Nothing(), use_reinit=true, - prefix="continuity_reinit", tspan=tspan) - -# Clip negative pressure to be able to use `SummationDensity` -state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure, - background_pressure=atmospheric_pressure, - clip_negative_pressure=true) + correction=nothing, use_reinit=true, + prefix="continuity_reinit", tspan=tspan, + fluid_density=fluid_density, density_diffusion=nothing) 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] + local smoothing_length = smoothing_length_dict[correction_name] println("="^100) println("fluid/dam_break_2d.jl with ", correction_name) @@ -48,9 +83,12 @@ for correction_name in keys(correction_dict) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length, - boundary_density_calculator=boundary_density_calculator, + boundary_density_calculator=SummationDensity(), fluid_density_calculator=fluid_density_calculator, correction=correction, use_reinit=false, - state_equation=state_equation, - prefix="$(correction_name)", tspan=tspan) + clip_negative_pressure=(fluid_density_calculator isa SummationDensity), + smoothing_kernel=smoothing_kernel, + prefix="$(correction_name)", tspan=tspan, + fluid_density=fluid_density, density_diffusion=Nothing(), + boundary_layers=5) end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c9bff5557..87789eeb1 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -5,7 +5,7 @@ using Reexport: @reexport using Dates using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect using FastPow: @fastpow -using LinearAlgebra: norm, dot, I, tr +using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using Morton: cartesian2morton using MuladdMacro: @muladd using Polyester: Polyester, @batch @@ -55,7 +55,8 @@ export examples_dir, trixi_include export trixi2vtk export RectangularTank, RectangularShape, SphereShape export VoxelSphere, RoundSphere, reset_wall! -export ShepardKernelCorrection, KernelGradientCorrection, AkinciFreeSurfaceCorrection +export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, + GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection export nparticles end # module diff --git a/src/general/corrections.jl b/src/general/corrections.jl index 595a61f0e..21fc8f89d 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -1,5 +1,4 @@ # Sorted in order of computational cost - @doc raw""" AkinciFreeSurfaceCorrection(rho0) @@ -11,6 +10,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 an SPH particle in the middle of a volume +at rest, its density will be identical to the rest density ``\rho_0``. If we now consider an SPH +particle at a free surface at rest, it will have neighbors missing in the direction normal to +the surface, which will result in a lower density. If we calculate the correction factor +```math +k = \rho_0/\rho_\text{mean}, +``` +this value will be about ~1.5 for particles at the free surface and can then be used to increase +the pressure and viscosity accordingly. + # Arguments - `rho0`: Rest density. @@ -71,19 +80,19 @@ to an improvement, especially at free surfaces. [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. + 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) - Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". - In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159-193. + In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159--193. [doi:10.1016/S0045-7825(96)01082-1](https://doi.org/10.1016/S0045-7825(96)01082-1) """ struct ShepardKernelCorrection end @doc raw""" - KernelGradientCorrection() + KernelCorrection() -Kernel gradient correction uses Shepard interpolation to obtain a 0-th order accurate result, which +Kernel correction uses Shepard interpolation to obtain a 0-th order accurate result, 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. @@ -93,7 +102,7 @@ c(x) = \sum_{b=1} V_b W_b(x) ``` The gradient of corrected kernel is determined by ```math -\nabla \tilde{W}_{b}(r) =\frac{\nabla W_{b}(r) - \gamma(r)}{\sum_{b=1} V_b W_b(r)} , \quad \text{where} \quad +\nabla \tilde{W}_{b}(r) =\frac{\nabla W_{b}(r) - W_b(r) \gamma(r)}{\sum_{b=1} V_b W_b(r)} , \quad \text{where} \quad \gamma(r) = \frac{\sum_{b=1} V_b \nabla W_b(r)}{\sum_{b=1} V_b W_b(r)}. ``` @@ -113,30 +122,64 @@ This correction can be applied with [`SummationDensity`](@ref) and [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. + 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) - Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159-193. [doi:10.1016/S0045-7825(96)01082-1](https://doi.org/10.1016/S0045-7825(96)01082-1) """ -struct KernelGradientCorrection end +struct KernelCorrection end -function kernel_correction_coefficient(system, particle) +@doc raw""" + MixedKernelGradientCorrection() + +Combines [`GradientCorrection`](@ref) and [`KernelCorrection`](@ref), +which results in a 1st-order-accurate SPH method. + +# Notes: +- Stability issues, 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 end + +function kernel_correction_coefficient(system::FluidSystem, particle) return system.cache.kernel_correction_coefficient[particle] end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - density_calculator, correction) +function kernel_correction_coefficient(system::BoundarySystem, particle) + return system.boundary_model.cache.kernel_correction_coefficient[particle] +end + +function compute_correction_values!(system, correction, v, u, v_ode, u_ode, semi, + density_calculator) return system end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - ::SummationDensity, ::ShepardKernelCorrection) +function compute_correction_values!(system, ::ShepardKernelCorrection, v, u, v_ode, u_ode, + semi, + ::SummationDensity) return compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, system.cache.kernel_correction_coefficient) end +function compute_correction_values!(system::BoundarySystem, ::ShepardKernelCorrection, v, u, + v_ode, u_ode, semi, + ::SummationDensity) + return compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, + system.boundary_model.cache.kernel_correction_coefficient) +end + function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, kernel_correction_coefficient) set_zero!(kernel_correction_coefficient) @@ -153,8 +196,9 @@ function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, # Loop over all pairs of particles and neighbors within the kernel cutoff for_particle_neighbor(system, neighbor_system, system_coords, - neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance + neighbor_coords, neighborhood_search, + particles=eachparticle(system)) do particle, neighbor, + pos_diff, distance rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) m_b = hydrodynamic_mass(neighbor_system, neighbor) volume = m_b / rho_b @@ -167,16 +211,43 @@ function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi, return kernel_correction_coefficient end -function dw_gamma(system, particle) +function dw_gamma(system::FluidSystem, particle) return extract_svector(system.cache.dw_gamma, system, particle) end -function compute_correction_values!(system, v, u, v_ode, u_ode, semi, - ::Union{SummationDensity, ContinuityDensity}, - ::KernelGradientCorrection) - (; cache) = system - (; kernel_correction_coefficient, dw_gamma) = cache +function dw_gamma(system::BoundarySystem, particle) + return extract_svector(system.boundary_model.cache.dw_gamma, system, particle) +end + +function compute_correction_values!(system::FluidSystem, + correction::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, + v_ode, u_ode, semi, + density_calculator) + compute_correction_values!(system, correction, + v, u, v_ode, u_ode, semi, + density_calculator, + system.cache.kernel_correction_coefficient, + system.cache.dw_gamma) +end +function compute_correction_values!(system::BoundarySystem, + correction::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, + v_ode, u_ode, semi, + density_calculator) + compute_correction_values!(system, correction, v, u, v_ode, u_ode, semi, + density_calculator, + system.boundary_model.cache.kernel_correction_coefficient, + system.boundary_model.cache.dw_gamma) +end + +function compute_correction_values!(system, + ::Union{KernelCorrection, + MixedKernelGradientCorrection}, v, u, v_ode, + u_ode, semi, + density_calculator, kernel_correction_coefficient, + dw_gamma) set_zero!(kernel_correction_coefficient) set_zero!(dw_gamma) @@ -193,7 +264,9 @@ function compute_correction_values!(system, v, u, v_ode, u_ode, semi, # Loop over all pairs of particles and neighbors within the kernel cutoff for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance + neighborhood_search, + particles=eachparticle(system)) do particle, neighbor, + pos_diff, distance rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) m_b = hydrodynamic_mass(neighbor_system, neighbor) volume = m_b / rho_b @@ -250,21 +323,46 @@ the correction matrix $\bm{L}_a$ is evaluated explicitly as !!! note - Stability issues arise, especially when particles separate into small clusters. - Doubles the computational effort. +- Better stability with smoother smoothing Kernels with larger support, e.g. [`SchoenbergQuinticSplineKernel`](@ref) or [`WendlandC6Kernel`](@ref). +- Set `dt_max =< 1e-3` for stability. ## 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. + 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. + 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 -function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, system, - coordinates, density_fun) +@doc raw""" + BlendedGradientCorrection() + +Calculate a blended gradient to reduce the stability issues of the [`GradientCorrection`](@ref). + +This calculates the following, +```math +\tilde\nabla A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i +``` +with ``0 \leq \lambda \leq 1`` being the blending factor. + +# Arguments +- `blending_factor`: Blending factor between corrected and regular SPH gradient. +""" +struct BlendedGradientCorrection{ELTYPE <: Real} + blending_factor::ELTYPE + + function BlendedGradientCorrection(blending_factor) + return new{eltype(blending_factor)}(blending_factor) + end +end + +# Called only by DensityDiffusion and TLSPH +function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, + system, coordinates, density_fun) (; mass) = system set_zero!(corr_matrix) @@ -275,12 +373,12 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s neighborhood_search; particles=eachparticle(system)) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0. - distance < sqrt(eps()) && return - volume = mass[neighbor] / density_fun(neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + + iszero(grad_kernel) && return + result = volume * grad_kernel * pos_diff' @inbounds for j in 1:ndims(system), i in 1:ndims(system) @@ -288,20 +386,95 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, s end end - @threaded for particle in eachparticle(system) - L = extract_smatrix(corr_matrix, system, particle) - result = inv(L) + correction_matrix_inversion_step!(corr_matrix, system) - 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 + return corr_matrix +end - @inbounds for j in 1:ndims(system), i in 1:ndims(system) - corr_matrix[i, j, particle] = result[i, j] +function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, + coordinates, v_ode, u_ode, semi, + correction, smoothing_length, smoothing_kernel) + set_zero!(corr_matrix) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + @trixi_timeit timer() "compute correction matrix" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + neighborhood_search = neighborhood_searches(system, neighbor_system, semi) + + for_particle_neighbor(system, neighbor_system, coordinates, neighbor_coords, + neighborhood_search; + particles=eachparticle(system)) do particle, + neighbor, + pos_diff, + distance + volume = hydrodynamic_mass(neighbor_system, neighbor) / + particle_density(v_neighbor_system, neighbor_system, neighbor) + + function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, + smoothing_length, system, particle) + return smoothing_kernel_grad(system, pos_diff, distance) + end + + # Compute gradient of corrected kernel + function compute_grad_kernel(correction::MixedKernelGradientCorrection, + smoothing_kernel, pos_diff, distance, + smoothing_length, system, particle) + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length, KernelCorrection(), system, + particle) + end + + grad_kernel = compute_grad_kernel(correction, smoothing_kernel, pos_diff, + distance, smoothing_length, system, particle) + + iszero(grad_kernel) && return + + L = volume * grad_kernel * pos_diff' + + # pos_diff is always x_a - x_b hence * -1 to switch the order to x_b - x_a + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] -= L[i, j] + end end end + correction_matrix_inversion_step!(corr_matrix, system) + return corr_matrix end + +function correction_matrix_inversion_step!(corr_matrix, system) + @threaded for particle in eachparticle(system) + L = extract_smatrix(corr_matrix, system, particle) + norm_ = norm(L) + + # The norm value is quasi-zero, so there are probably no neighbors for this particle + if norm_ < sqrt(eps()) + # The correction matrix is set to an identity matrix, which effectively disables + # the correction for this particle. + @inbounds for j in 1:ndims(system), i in 1:ndims(system) + corr_matrix[i, j, particle] = 0.0 + end + @inbounds for i in 1:ndims(system) + corr_matrix[i, i, particle] = 1.0 + end + continue + end + + det_ = abs(det(L)) + @fastmath if det_ < 1e-6 * norm_ + 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 + else + 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 +end diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 4383d55e3..732a5eb2a 100644 --- a/src/general/density_calculators.jl +++ b/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, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) diff --git a/src/general/general.jl b/src/general/general.jl index e5cf8bdd4..a5c3f8fc5 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -1,3 +1,13 @@ +abstract type System{NDIMS} end + +abstract type FluidSystem{NDIMS} <: System{NDIMS} end + +abstract type SolidSystem{NDIMS} <: System{NDIMS} end + +abstract type BoundarySystem{NDIMS} <: System{NDIMS} end + +timer_name(::FluidSystem) = "fluid" + @inline function set_zero!(du) du .= zero(eltype(du)) diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index c71869ee5..c2f48ebba 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -52,6 +52,8 @@ abstract type SmoothingKernel{NDIMS} end @inline Base.ndims(::SmoothingKernel{NDIMS}) where {NDIMS} = NDIMS @inline function kernel_grad(kernel, pos_diff, distance, h) + distance < sqrt(eps()) && return zero(pos_diff) + return kernel_deriv(kernel, distance, h) / distance * pos_diff end @@ -60,12 +62,37 @@ end return kernel_grad(kernel, pos_diff, distance, h) end -@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, - ::KernelGradientCorrection, system, particle) - return (kernel_grad(kernel, pos_diff, distance, h) .- dw_gamma(system, particle)) / +@inline function corrected_kernel_grad(kernel_, pos_diff, distance, h, ::KernelCorrection, + system, particle) + return (kernel_grad(kernel_, pos_diff, distance, h) .- + kernel(kernel_, distance, h) * dw_gamma(system, particle)) / kernel_correction_coefficient(system, particle) end +@inline function corrected_kernel_grad(kernel, pos_diff, distance, h, + corr::BlendedGradientCorrection, system, + particle) + (; blending_factor) = corr + + grad = kernel_grad(kernel, pos_diff, distance, h) + return (1 - blending_factor) * grad + + blending_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, KernelCorrection(), + system, particle) + return correction_matrix(system, particle) * grad +end + @doc raw""" GaussianKernel{NDIMS}() diff --git a/src/general/system.jl b/src/general/system.jl index a65796a66..b27b229e3 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -1,5 +1,3 @@ -abstract type System{NDIMS} end - initialize!(system, neighborhood_search) = system @inline Base.ndims(::System{NDIMS}) where {NDIMS} = NDIMS @@ -80,6 +78,12 @@ end return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) end +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance) + (; smoothing_kernel, smoothing_length) = system.boundary_model + + return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) +end + # This is dispatched for some system types @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 73432d9c7..4ca85c33e 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -74,7 +74,7 @@ We provide five options to compute the boundary density and pressure, determined In: Computers, Materials and Continua 5 (2007), pages 173-184. [doi: 10.3970/cmc.2007.005.173](https://doi.org/10.3970/cmc.2007.005.173) """ -struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, C} +struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, COR, C} pressure :: Vector{ELTYPE} hydrodynamic_mass :: Vector{ELTYPE} state_equation :: SE @@ -82,37 +82,34 @@ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, SE, K, V, C} smoothing_kernel :: K smoothing_length :: ELTYPE viscosity :: V + correction :: COR cache :: C function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=NoViscosity(), - state_equation=nothing) + state_equation=nothing, correction=nothing) pressure = initial_boundary_pressure(initial_density, density_calculator, state_equation) + NDIMS = ndims(smoothing_kernel) n_particles = length(initial_density) - cache = (; create_cache_model(viscosity, n_particles, ndims(smoothing_kernel))..., + cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., create_cache_model(initial_density, density_calculator)...) + cache = (; + create_cache_model(correction, initial_density, NDIMS, + n_particles)..., cache...) new{typeof(density_calculator), eltype(initial_density), typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), - typeof(cache)}(pressure, hydrodynamic_mass, state_equation, density_calculator, - smoothing_kernel, smoothing_length, viscosity, cache) + typeof(correction), typeof(cache)}(pressure, hydrodynamic_mass, state_equation, + density_calculator, + smoothing_kernel, smoothing_length, + viscosity, correction, cache) end end -function Base.show(io::IO, model::BoundaryModelDummyParticles) - @nospecialize model # reduce precompilation time - - print(io, "BoundaryModelDummyParticles(") - print(io, model.density_calculator |> typeof |> nameof) - print(io, ", ") - print(io, model.viscosity |> typeof |> nameof) - print(io, ")") -end - @doc raw""" AdamiPressureExtrapolation() @@ -188,39 +185,27 @@ reference pressure (the corresponding pressure to the reference density by the s """ struct PressureZeroing end -# For most density calculators, the pressure is updated in every step -initial_boundary_pressure(initial_density, density_calculator, _) = similar(initial_density) -# Pressure mirroring does not use the pressure, so we set it to zero for the visualization -initial_boundary_pressure(initial_density, ::PressureMirroring, _) = zero(initial_density) +create_cache_model(correction, density, NDIMS, nparticles) = (;) -# For pressure zeroing, set the pressure to the reference pressure (zero with free surfaces) -function initial_boundary_pressure(initial_density, ::PressureZeroing, state_equation) - return state_equation.(initial_density) +function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density)) end -# With EDAC, just use zero pressure -function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing) - return zero(initial_density) +function create_cache_model(::KernelCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density), dw_gamma) end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, smoothing_length, - grad_kernel, - boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, - fluid_density_calculator) - - # Use `p_a` as pressure for both particles with `PressureMirroring` - return pressure_acceleration(pressure_correction, m_b, p_a, p_a, rho_a, rho_b, - grad_kernel, fluid_density_calculator) +function create_cache_model(::Union{GradientCorrection, BlendedGradientCorrection}, density, + NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; correction_matrix) end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, smoothing_length, - grad_kernel, - boundary_model::BoundaryModelDummyParticles, - fluid_density_calculator) - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, fluid_density_calculator) +function create_cache_model(::MixedKernelGradientCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end function create_cache_model(initial_density, @@ -264,6 +249,90 @@ function reset_cache!(cache, viscosity::ViscosityAdami) return cache end +function Base.show(io::IO, model::BoundaryModelDummyParticles) + @nospecialize model # reduce precompilation time + + print(io, "BoundaryModelDummyParticles(") + print(io, model.density_calculator |> typeof |> nameof) + print(io, ", ") + print(io, model.viscosity |> typeof |> nameof) + print(io, ")") +end + +# For most density calculators, the pressure is updated in every step +initial_boundary_pressure(initial_density, density_calculator, _) = similar(initial_density) +# Pressure mirroring does not use the pressure, so we set it to zero for the visualization +initial_boundary_pressure(initial_density, ::PressureMirroring, _) = zero(initial_density) + +# For pressure zeroing, set the pressure to the reference pressure (zero with free surfaces) +function initial_boundary_pressure(initial_density, ::PressureZeroing, state_equation) + return state_equation.(initial_density) +end + +# With EDAC, just use zero pressure +function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing) + return zero(initial_density) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, + fluid_density_calculator, + correction) + + # Use `p_a` as pressure for both particles with `PressureMirroring` + return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_a, rho_a, rho_b, + grad_kernel, fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles, + fluid_density_calculator, + correction) + return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, + grad_kernel, fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, W_a, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles{<:PressureMirroring}, + fluid_density_calculator, + correction::Union{KernelCorrection, + 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, W_a, W_b, + fluid_density_calculator) +end + +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, W_a, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelDummyParticles, + fluid_density_calculator, + correction::Union{KernelCorrection, + 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, W_a, W_b, + fluid_density_calculator) +end + @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle) return particle_density(v, model.density_calculator, model, particle) end @@ -302,15 +371,55 @@ end @inline function update_pressure!(boundary_model::BoundaryModelDummyParticles, system, v, u, v_ode, u_ode, semi) - (; density_calculator) = boundary_model + (; density_calculator, correction) = boundary_model + + compute_correction_values!(system, + correction, v, u, v_ode, u_ode, semi, density_calculator) + + 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, + density_calculator) compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) return boundary_model end -function compute_density!(boundary_model, ::SummationDensity, - system, v, u, v_ode, u_ode, semi) +function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, + correction, density_calculator) + return boundary_model +end + +function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, + corr::ShepardKernelCorrection, ::SummationDensity) + boundary_model.cache.density ./= boundary_model.cache.kernel_correction_coefficient +end + +function compute_gradient_correction_matrix!(correction, boundary_model, system, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + boundary_model, + 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, + v_ode, u_ode, semi, correction, smoothing_length, + smoothing_kernel) +end + +function compute_density!(boundary_model, ::SummationDensity, system, v, u, v_ode, u_ode, + semi) (; cache) = boundary_model (; density) = cache # Density is in the cache for SummationDensity @@ -487,3 +596,15 @@ end # The density is constant when using EDAC return density end + +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, + distance, particle) + (; smoothing_kernel, smoothing_length, correction) = system.boundary_model + + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length, correction, system, particle) +end + +@inline function correction_matrix(system::BoundarySystem, particle) + extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) +end diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index e608128a6..fa476e530 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -76,11 +76,13 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff::SVector{NDIMS}, - smoothing_length, grad_kernel, - boundary_model::BoundaryModelMonaghanKajtar, - density_calculator) where {NDIMS} +@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff::SVector{NDIMS}, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelMonaghanKajtar, + fluid_density_calculator, + correction) where {NDIMS} (; K, beta, boundary_particle_spacing) = boundary_model distance = norm(pos_diff) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index e3cc8c29d..a5be4ab93 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -7,7 +7,7 @@ The interaction between fluid and boundary particles is specified by the boundar For moving boundaries, a [`BoundaryMovement`](@ref) can be passed with the keyword argument `movement`. """ -struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: System{NDIMS} +struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: BoundarySystem{NDIMS} coordinates :: Array{ELTYPE, 2} boundary_model :: BM movement :: M @@ -266,13 +266,17 @@ function viscosity_model(system::BoundarySPHSystem) end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system::BoundarySPHSystem, - density_calculator) + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, + neighbor_system::BoundarySPHSystem, + density_calculator, correction) (; boundary_model) = neighbor_system (; smoothing_length) = particle_system - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - smoothing_length, grad_kernel, boundary_model, - density_calculator) + return pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model, + density_calculator, correction) end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 844e0a91d..712a6c7dc 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -1,7 +1,3 @@ -abstract type FluidSystem{NDIMS} <: System{NDIMS} end - -timer_name(::FluidSystem) = "fluid" - @inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] function write_u0!(u0, system::FluidSystem) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 11a9e8d4e..f0461f477 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -21,9 +21,6 @@ function interact!(dv, v_particle_system, u_particle_system, for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return - rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = 0.5 * (rho_a + rho_b) @@ -42,9 +39,9 @@ function interact!(dv, v_particle_system, u_particle_system, p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) dv_pressure = pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system, - density_calculator) + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, neighbor_system, + density_calculator, correction) dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, @@ -58,6 +55,7 @@ function interact!(dv, v_particle_system, u_particle_system, # debug_array[i, particle] += dv_pressure[i] end + # TODO If variable smoothing_length is used, this should use the neighbor smoothing length continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, particle_system, neighbor_system, grad_kernel) @@ -73,13 +71,33 @@ function interact!(dv, v_particle_system, u_particle_system, end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, neighbor_system, - density_calculator) + rho_a, rho_b, pos_diff, distance, + grad_kernel, particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + density_calculator, + correction) + + # 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 - # By default, just call the pressure acceleration formulation corresponding to the density calculator - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, density_calculator) +@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + W_a, particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + density_calculator, + correction::Union{KernelCorrection, + GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}) + W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor) + + # 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, W_a, W_b, density_calculator) end # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic @@ -87,8 +105,10 @@ end # 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. -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, ::ContinuityDensity) +# 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) return (-m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel) * pressure_correction end @@ -97,11 +117,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. -@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, - grad_kernel, ::SummationDensity) +# 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) 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, W_a, 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, W_a, 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! @inline function continuity_equation!(dv, density_calculator::SummationDensity, v_particle_system, v_neighbor_system, @@ -111,6 +147,7 @@ end return dv end +# This formulation was chosen to be consistent with the used pressure_acceleration formulations. @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -136,6 +173,9 @@ end particle_system::WeaklyCompressibleSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem, grad_kernel) + # Density diffusion terms are all zero for distance zero + distance < sqrt(eps()) && return + (; delta) = density_diffusion (; smoothing_length, state_equation) = particle_system (; sound_speed) = state_equation diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index e65d1465a..4108466a6 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -87,12 +87,25 @@ end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) function create_cache_wcsph(::ShepardKernelCorrection, density, NDIMS, n_particles) - (; kernel_correction_coefficient=similar(density)) + return (; kernel_correction_coefficient=similar(density)) end -function create_cache_wcsph(::KernelGradientCorrection, density, NDIMS, n_particles) +function create_cache_wcsph(::KernelCorrection, 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_wcsph(::Union{GradientCorrection, BlendedGradientCorrection}, density, + NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; correction_matrix) +end + +function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + + return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end function create_cache_wcsph(n_particles, ELTYPE, ::SummationDensity) @@ -162,7 +175,7 @@ initialize!(system::WeaklyCompressibleSPHSystem, neighborhood_search) = system function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, density_diffusion) = system + (; density_calculator, density_diffusion, correction) = system compute_density!(system, u, u_ode, semi, density_calculator) @@ -188,8 +201,11 @@ end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, correction) = system - compute_correction_values!(system, v, u, v_ode, u_ode, semi, - density_calculator, correction) + compute_correction_values!(system, + correction, v, u, v_ode, u_ode, semi, density_calculator) + + 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, density_calculator) @@ -198,17 +214,37 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end -function kernel_correct_density!(system, v, u, v_ode, u_ode, semi, - correction, density_calculator) +function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, + semi, correction, density_calculator) return system end -function kernel_correct_density!(system, v, u, v_ode, u_ode, semi, - ::Union{ShepardKernelCorrection, KernelGradientCorrection}, - ::SummationDensity) +function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, + semi, corr::ShepardKernelCorrection, ::SummationDensity) system.cache.density ./= system.cache.kernel_correction_coefficient end +function compute_gradient_correction_matrix!(correction, + system::WeaklyCompressibleSPHSystem, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + system::WeaklyCompressibleSPHSystem, u, + 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, + v_ode, u_ode, semi, correction, smoothing_length, + smoothing_kernel) +end + function reinit_density!(vu_ode, semi) v_ode, u_ode = vu_ode.x @@ -313,3 +349,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 diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index bd9505505..f5c17dd68 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -61,7 +61,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem) (; boundary_model) = particle_system - (; density_calculator, state_equation, viscosity) = neighbor_system + (; density_calculator, state_equation, viscosity, correction) = neighbor_system (; sound_speed) = state_equation system_coords = current_coordinates(u_particle_system, particle_system) @@ -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,10 @@ 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, rho_b, rho_a, pos_diff, - neighbor_system.smoothing_length, - grad_kernel, boundary_model, density_calculator) + dv_boundary = pressure_acceleration(1.0, m_a, p_b, p_a, + rho_b, rho_a, pos_diff, distance, grad_kernel, + neighbor_system, neighbor, particle_system, + density_calculator, correction) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 013048332..14a4b37ab 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -75,7 +75,7 @@ The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGa In: International Journal for Numerical Methods in Engineering 48 (2000), pages 1359–1400. [doi: 10.1002/1097-0207](https://doi.org/10.1002/1097-0207) """ -struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIMS} +struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} initial_coordinates :: Array{ELTYPE, 2} # [dimension, particle] current_coordinates :: Array{ELTYPE, 2} # [dimension, particle] @@ -94,7 +94,6 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM acceleration :: SVector{NDIMS, ELTYPE} boundary_model :: BM penalty_force :: PF - function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio, boundary_model; @@ -130,17 +129,15 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM ((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, + 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) + lame_lambda, lame_mu, smoothing_kernel, + smoothing_length, acceleration_, boundary_model, + penalty_force) end end @@ -234,6 +231,7 @@ end @inline function correction_matrix(system, particle) extract_smatrix(system.correction_matrix, system, particle) end + @inline function deformation_gradient(system, particle) extract_smatrix(system.deformation_grad, system, particle) end @@ -412,16 +410,18 @@ function viscosity_model(system::TotalLagrangianSPHSystem) end @inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b, - rho_a, rho_b, pos_diff, grad_kernel, - particle_system, + rho_a, rho_b, pos_diff, distance, grad_kernel, + particle_system, neighbor, neighbor_system::TotalLagrangianSPHSystem, - density_calculator) + density_calculator, correction) (; boundary_model) = neighbor_system (; smoothing_length) = particle_system # Pressure acceleration for fluid-solid interaction. This is identical to # `pressure_acceleration` for the `BoundarySPHSystem`. - return pressure_acceleration(pressure_correction, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - smoothing_length, grad_kernel, boundary_model, - density_calculator) + return pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, distance, + smoothing_length, grad_kernel, + particle_system, neighbor, neighbor_system, + boundary_model, density_calculator, correction) end diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index e7ee22392..18f98f117 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -6,14 +6,14 @@ integral_2d_radial, _ = quadgk(r -> r * TrixiParticles.kernel(smk, r, 1.0), 0, TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) - return 2 * π * integral_2d_radial + return 2 * pi * integral_2d_radial end function integrate_kernel_3d(smk) integral_3d_radial, _ = quadgk(r -> r^2 * TrixiParticles.kernel(smk, r, 1.0), 0, TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) - return 4 * π * integral_3d_radial + return 4 * pi * integral_3d_radial end # Treat the truncated Gaussian kernel separately diff --git a/test/schemes/fluid/weakly_compressible_sph/rhs.jl b/test/schemes/fluid/weakly_compressible_sph/rhs.jl index 6e98f1c2f..a4a6a5e8a 100644 --- a/test/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -51,13 +51,17 @@ # Compute accelerations a -> b and b -> a dv1 = TrixiParticles.pressure_acceleration(1.0, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - grad_kernel, system, system, - density_calculator) + distance, + grad_kernel, system, -1, + system, + density_calculator, nothing) dv2 = TrixiParticles.pressure_acceleration(1.0, m_a, p_b, p_a, rho_b, rho_a, -pos_diff, - -grad_kernel, system, system, - density_calculator) + distance, + -grad_kernel, system, -1, + system, + density_calculator, nothing) # Test that both forces are identical but in opposite directions @test isapprox(m_a * dv1, -m_b * dv2, rtol=2eps()) diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 99d0795bc..98ab834b6 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -88,7 +88,7 @@ "SphereShape 2D", "RectangularShape 2D with ShepardKernelCorrection", "RectangularShape 2D with AkinciFreeSurfaceCorrection", - "RectangularShape 2D with KernelGradientCorrection", + "RectangularShape 2D with KernelCorrection", ] NDIMS_ = [2, 3, 2, 3, 2, 2, 2, 2] density_calculators = [ @@ -103,7 +103,7 @@ Nothing(), ShepardKernelCorrection(), AkinciFreeSurfaceCorrection(1000.0), - KernelGradientCorrection(), + KernelCorrection(), ] @testset "$(setup_names[i])" for i in eachindex(setups) @@ -146,7 +146,7 @@ if density_calculator isa SummationDensity @test length(system.cache.density) == size(setup.coordinates, 2) end - if corr isa ShepardKernelCorrection || corr isa KernelGradientCorrection + if corr isa ShepardKernelCorrection || corr isa KernelCorrection @test length(system.cache.kernel_correction_coefficient) == size(setup.coordinates, 2) end