diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 366cb1183..a780e9751 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.16.23 + uses: crate-ci/typos@v1.16.26 diff --git a/.gitignore b/.gitignore index b83d9967c..e017a301c 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,10 @@ docs/src/license.md docs/src/code_of_conduct.md docs/src/news.md run +*.json +*.png +*.gif +*.svg +*.vtu +*.pvd +*.mp4 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/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index cb9a8d05b..6dbdf5927 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -101,7 +101,14 @@ function symplectic_euler!(velocities, coordinates, semi) end # One RHS evaluation is so fast that timers make it multiple times slower. -TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) +# Disabling timers throws a warning, which we suppress here in order to make the tests pass. +# For some reason, this only works with a file and not with devnull. See issue #332. +filename = tempname() +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + end +end @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi)) @@ -112,5 +119,9 @@ end @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi)) -# Enable timers again. -TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) +# Enable timers again +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + end +end diff --git a/examples/n_body/n_body_solar_system.jl b/examples/n_body/n_body_solar_system.jl index f063d0599..a14114fb0 100644 --- a/examples/n_body/n_body_solar_system.jl +++ b/examples/n_body/n_body_solar_system.jl @@ -44,7 +44,14 @@ saving_callback = SolutionSavingCallback(dt=10day) callbacks = CallbackSet(info_callback, saving_callback) # One RHS evaluation is so fast that timers make it multiple times slower. -TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) +# Disabling timers throws a warning, which we suppress here in order to make the tests pass. +# For some reason, this only works with a file and not with devnull. See issue #332. +filename = tempname() +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + end +end sol = solve(ode, SymplecticEuler(), dt=1.0e5, @@ -53,5 +60,9 @@ sol = solve(ode, SymplecticEuler(), @printf("%.9e\n", energy(ode.u0.x..., particle_system, semi)) @printf("%.9e\n", energy(sol.u[end].x..., particle_system, semi)) -# Enable timers again. -TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) +# Enable timers again +open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + end +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/semidiscretization.jl b/src/general/semidiscretization.jl index 6c3bf4796..6cfe1f8d2 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -168,8 +168,7 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`. function semidiscretize(semi, tspan; reset_threads=true) (; systems) = semi - @assert all(system -> eltype(system) === eltype(systems[1]), - systems) + @assert all(system -> eltype(system) === eltype(systems[1]), systems) ELTYPE = eltype(systems[1]) # Optionally reset Polyester.jl threads. See @@ -190,10 +189,8 @@ function semidiscretize(semi, tspan; reset_threads=true) end end - sizes_u = (u_nvariables(system) * n_moving_particles(system) - for system in systems) - sizes_v = (v_nvariables(system) * n_moving_particles(system) - for system in systems) + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in systems) + sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in systems) u0_ode = Vector{ELTYPE}(undef, sum(sizes_u)) v0_ode = Vector{ELTYPE}(undef, sum(sizes_v)) 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 c967b827c..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 @@ -319,18 +428,24 @@ end function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityDensity}, system, v, u, v_ode, u_ode, semi) - (; state_equation, pressure) = boundary_model # Limit pressure to be non-negative to avoid attractive forces between fluid and # boundary particles at free surfaces (sticking artifacts). - @trixi_timeit timer() "state equation" @threaded for particle in eachparticle(system) - pressure[particle] = max(state_equation(particle_density(v, boundary_model, - particle)), 0.0) + @threaded for particle in eachparticle(system) + apply_state_equation!(boundary_model, particle_density(v, boundary_model, + particle), particle) end return boundary_model end +# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). +# Otherwise, `@threaded` does not work here with Julia ARM on macOS. +# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. +@inline function apply_state_equation!(boundary_model, density, particle) + boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) +end + function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, v_ode, u_ode, semi) (; pressure, state_equation, cache, viscosity) = boundary_model @@ -358,19 +473,29 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, end @trixi_timeit timer() "inverse state equation" @threaded for particle in eachparticle(system) - # The summation is only over fluid particles, thus the volume stays zero when a boundary - # particle isn't surrounded by fluid particles. - # Check the volume to avoid NaNs in pressure and velocity. - if volume[particle] > eps() - pressure[particle] /= volume[particle] - - # To impose no-slip condition - compute_wall_velocity!(viscosity, system, system_coords, particle) - end - - # Apply inverse state equation to compute density (not used with EDAC) - inverse_state_equation!(density, state_equation, pressure, particle) + compute_adami_density!(boundary_model, system, system_coords, particle) + end +end + +# Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). +# Otherwise, `@threaded` does not work here with Julia ARM on macOS. +# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. +function compute_adami_density!(boundary_model, system, system_coords, particle) + (; pressure, state_equation, cache, viscosity) = boundary_model + (; volume, density) = cache + + # The summation is only over fluid particles, thus the volume stays zero when a boundary + # particle isn't surrounded by fluid particles. + # Check the volume to avoid NaNs in pressure and velocity. + if volume[particle] > eps() + pressure[particle] /= volume[particle] + + # To impose no-slip condition + compute_wall_velocity!(viscosity, system, system_coords, particle) end + + # Apply inverse state equation to compute density (not used with EDAC) + inverse_state_equation!(density, state_equation, pressure, particle) end function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZeroing}, @@ -471,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 778999913..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 @@ -251,7 +287,8 @@ end # Use this function to avoid passing closures to Polyester.jl with `@batch` (`@threaded`). # Otherwise, `@threaded` does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. -@inline function apply_state_equation!(system, density, particle) +@inline function apply_state_equation!(system::WeaklyCompressibleSPHSystem, density, + particle) system.pressure[particle] = system.state_equation(density) end @@ -312,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/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index e38f5cdee..aaf013ada 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -1,6 +1,7 @@ @doc raw""" RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density; n_layers=1, spacing_ratio=1.0, + min_coordinates=zeros(length(fluid_size)), init_velocity=zeros(length(fluid_size)), boundary_density=fluid_density, faces=Tuple(trues(2 * length(fluid_size)))) @@ -18,6 +19,7 @@ Rectangular tank filled with a fluid to set up dam-break-style simulations. - `spacing_ratio`: Ratio of `particle_spacing` to boundary particle spacing. A value of 2 means that the boundary particle spacing will be half the fluid particle spacing. +- `min_coordinates`: Coordinates of the corner in negative coordinate directions. - `init_velocity`: The initial velocity of each fluid particle as `(x, y)` (or `(x, y, z)` in 3D). - `boundary_density`: Density of each boundary particle (by default set to the rest density) - `faces`: By default all faces are generated. Set faces by passing a @@ -58,6 +60,7 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} function RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density; pressure=0.0, n_layers=1, spacing_ratio=1.0, + min_coordinates=zeros(length(fluid_size)), init_velocity=zeros(length(fluid_size)), boundary_density=fluid_density, faces=Tuple(trues(2 * length(fluid_size))), @@ -119,6 +122,10 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} boundary_masses, boundary_densities, particle_spacing=boundary_spacing) + # Move the tank corner in the negative coordinate directions to the desired position + fluid.coordinates .+= min_coordinates + boundary.coordinates .+= min_coordinates + return new{NDIMS, 2 * NDIMS, ELTYPE}(fluid, boundary, fluid_size_, tank_size_, faces, face_indices, particle_spacing, spacing_ratio, n_layers, @@ -152,7 +159,7 @@ function fluid_particles_per_dimension(size::NTuple{3}, particle_spacing) "fluid length in x-direction") n_particles_y, new_y_size = round_n_particles(size[2], particle_spacing, "fluid length in y-direction") - n_particles_z, new_z_size = round_n_particles(size[2], particle_spacing, + n_particles_z, new_z_size = round_n_particles(size[3], particle_spacing, "fluid length in z-direction") return (n_particles_x, n_particles_y, n_particles_z), diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index 73a40a376..a2ab46b9c 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -127,15 +127,38 @@ as it will have corners (like a sphere in Minecraft). struct VoxelSphere end """ - RoundSphere() + RoundSphere(; start_angle=0.0, end_angle=2π) -Construct a sphere by nesting perfectly round concentric spheres. +Construct a sphere (or sphere segment) by nesting perfectly round concentric spheres. The resulting ball will be perfectly round, but will not have a regular inner structure. +# Keywords +- `start_angle`: The starting angle of the sphere segment in radians. It determines the + beginning point of the segment. The default is set to `0.0` representing + the positive x-axis. +- `end_angle`: The ending angle of the sphere segment in radians. It defines the termination + point of the segment. The default is set to `2pi`, completing a full sphere. + !!! note "Usage" See [`SphereShape`](@ref) on how to use this. + +!!! warning "Warning" + The sphere segment is intended for 2D geometries and hollow spheres. If used for filled + spheres or in a 3D context, results may not be accurate. """ -struct RoundSphere end +struct RoundSphere{AR} + angle_range::AR + + function RoundSphere(; start_angle=0.0, end_angle=2pi) + if start_angle > end_angle + throw(ArgumentError("`end_angle` should be greater than `start_angle`")) + end + + angle_range = (start_angle, end_angle) + + new{typeof(angle_range)}(angle_range) + end +end function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_position, n_layers, layer_outwards, tlsph) @@ -193,7 +216,7 @@ function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_pos return reinterpret(reshape, ELTYPE, coords) end -function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, +function sphere_shape_coords(sphere::RoundSphere, particle_spacing, radius, center, n_layers, layer_outwards, tlsph) if n_layers > 0 if layer_outwards @@ -227,7 +250,7 @@ function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, coords = zeros(length(center), 0) for layer in 0:(n_layers - 1) - sphere_coords = round_sphere(particle_spacing, + sphere_coords = round_sphere(sphere, particle_spacing, inner_radius + layer * particle_spacing, center) coords = hcat(coords, sphere_coords) end @@ -235,8 +258,11 @@ function sphere_shape_coords(::RoundSphere, particle_spacing, radius, center, return coords end -function round_sphere(particle_spacing, radius, center::SVector{2}) - n_particles = round(Int, 2pi * radius / particle_spacing) +function round_sphere(sphere, particle_spacing, radius, center::SVector{2}) + (; angle_range) = sphere + + theta = angle_range[2] - angle_range[1] + n_particles = round(Int, theta * radius / particle_spacing) if n_particles <= 2 # 2 or less particles produce weird, asymmetric results. @@ -244,8 +270,12 @@ function round_sphere(particle_spacing, radius, center::SVector{2}) return collect(reshape(center, (2, 1))) end - # Remove the last particle at 2pi, which overlaps with the first at 0 - t = LinRange(0, 2pi, n_particles + 1)[1:(end - 1)] + if !isapprox(theta, 2pi) + t = LinRange(angle_range[1], angle_range[2], n_particles + 1) + else + # Remove the last particle at 2pi, which overlaps with the first at 0 + t = LinRange(angle_range[1], angle_range[2], n_particles + 1)[1:(end - 1)] + end particle_coords = Array{Float64, 2}(undef, 2, length(t)) @@ -256,7 +286,7 @@ function round_sphere(particle_spacing, radius, center::SVector{2}) return particle_coords end -function round_sphere(particle_spacing, radius, center::SVector{3}) +function round_sphere(sphere, particle_spacing, radius, center::SVector{3}) # The number of particles can either be calculated in 2D or in 3D. # Let δ be the particle spacing and r the sphere radius. # @@ -374,7 +404,7 @@ function round_sphere(particle_spacing, radius, center::SVector{3}) circle_spacing = 1.0 end - circle_coords_2d = round_sphere(circle_spacing, circle_radius, + circle_coords_2d = round_sphere(sphere, circle_spacing, circle_radius, SVector(center[1], center[2])) circle_coords_3d = vcat(circle_coords_2d, center[3] .+ z * ones(1, size(circle_coords_2d, 2))) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 2f78d9934..8b2227a0b 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -70,7 +70,8 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix collection_file = joinpath(output_directory, add_opt_str_pre(prefix) * "$system_name") - pvd = paraview_collection(collection_file; append=true) + # Reset the collection when the iteration is 0 + pvd = paraview_collection(collection_file; append=iter > 0) points = periodic_coords(current_coordinates(u, system), periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] diff --git a/test/Project.toml b/test/Project.toml index 85ced8502..897dd26f0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] 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..09819e898 100644 --- a/test/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -38,9 +38,8 @@ TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 smoothing_length = -1.0 - initial_condition = InitialCondition(coordinates, velocity, mass, - density) - system = WeaklyCompressibleSPHSystem(initial_condition, + fluid = InitialCondition(coordinates, velocity, mass, density) + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, state_equation, smoothing_kernel, smoothing_length) @@ -51,13 +50,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()) @@ -66,142 +69,120 @@ end end - @testset verbose=true "`interact!`" begin - # The following tests for linear and angular momentum and total energy conservation - # are based on Sections 3.3.4 and 3.4.2 of - # Daniel J. Price. "Smoothed Particle Hydrodynamics and Magnetohydrodynamics." - # In: Journal of Computational Physics 231.3 (2012), pages 759–94. - # https://doi.org/10.1016/j.jcp.2010.12.011 - @testset verbose=true "Momentum and Total Energy Conservation" begin - # We are testing the momentum conservation of SPH with random initial configurations - density_calculators = [ContinuityDensity(), SummationDensity()] - - # Random initial configuration - mass = [ - [3.11, 1.55, 2.22, 3.48, 0.21, 3.73, 0.21, 3.45], - [0.82, 1.64, 1.91, 0.02, 0.08, 1.58, 4.94, 0.7], - ] - density = [ - [914.34, 398.36, 710.22, 252.54, 843.81, 694.73, 670.5, 539.14], - [280.15, 172.25, 267.1, 130.42, 382.3, 477.21, 848.31, 188.62], - ] - pressure = [ - [91438.0, 16984.0, 58638.0, 10590.0, 92087.0, 66586.0, 64723.0, 49862.0], - [31652.0, -21956.0, 2874.0, -12489.0, 27206.0, 32225.0, 42848.0, 3001.0], - ] - coordinates = [ - [0.16 0.55 0.08 0.58 0.52 0.26 0.32 0.99; - 0.76 0.6 0.47 0.4 0.25 0.79 0.45 0.63], - [0.4 0.84 0.47 0.02 0.64 0.85 0.02 0.15; - 0.83 0.62 0.99 0.57 0.25 0.72 0.34 0.69], - ] - velocity = [ - [1.05 0.72 0.12 1.22 0.67 0.85 1.42 -0.57; - 1.08 0.68 0.74 -0.27 -1.22 0.43 1.41 1.25], - [-1.84 -1.4 5.21 -5.99 -5.02 9.5 -4.51 -8.28; - 0.78 0.1 9.67 8.46 9.29 5.18 -4.83 -4.87]] - - # The state equation is only needed to unpack `sound_speed`, so we can mock - # it by using a `NamedTuple`. - state_equation = (; sound_speed=0.0) - smoothing_kernel = SchoenbergCubicSplineKernel{2}() - smoothing_length = 0.3 - search_radius = TrixiParticles.compact_support(smoothing_kernel, - smoothing_length) - - @testset "`$(nameof(typeof(density_calculator)))`" for density_calculator in density_calculators - for i in eachindex(mass) - initial_condition = InitialCondition(coordinates[i], velocity[i], - mass[i], density[i]) - system = WeaklyCompressibleSPHSystem(initial_condition, - density_calculator, - state_equation, smoothing_kernel, - smoothing_length) - - # Overwrite `system.pressure` - system.pressure .= pressure[i] - - u = coordinates[i] - if density_calculator isa SummationDensity - # Density is stored in the cache - v = velocity[i] - system.cache.density .= density[i] - else - # Density is integrated with `ContinuityDensity` - v = vcat(velocity[i], density[i]') - end + # The following tests for linear and angular momentum and total energy conservation + # are based on Sections 3.3.4 and 3.4.2 of + # Daniel J. Price. "Smoothed Particle Hydrodynamics and Magnetohydrodynamics." + # In: Journal of Computational Physics 231.3 (2012), pages 759–94. + # https://doi.org/10.1016/j.jcp.2010.12.011 + @testset verbose=true "Momentum and Total Energy Conservation" begin + # We are testing the momentum conservation of SPH with random initial configurations + density_calculators = [ContinuityDensity(), SummationDensity()] + + particle_spacing = 0.1 + + # The state equation is only needed to unpack `sound_speed`, so we can mock + # it by using a `NamedTuple`. + state_equation = (; sound_speed=0.0) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.2particle_spacing + search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length) + + @testset "`$(nameof(typeof(density_calculator)))`" for density_calculator in density_calculators + # Run three times with different seed for the random initial condition + for seed in 1:3 + # A larger number of particles will increase accumulated errors in the + # summation. A larger tolerance has to be used for the tests below. + fluid = rectangular_patch(particle_spacing, (3, 3), seed=seed) + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, + state_equation, smoothing_kernel, + smoothing_length) + n_particles = TrixiParticles.nparticles(system) + n_vars = TrixiParticles.v_nvariables(system) + + # Overwrite `system.pressure` because we skip the update step + system.pressure .= fluid.pressure + + u = fluid.coordinates + if density_calculator isa SummationDensity + # Density is stored in the cache + v = fluid.velocity + system.cache.density .= fluid.density + else + # Density is integrated with `ContinuityDensity` + v = vcat(fluid.velocity, fluid.density') + end - nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius, - TrixiParticles.eachparticle(system)) + nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius, + TrixiParticles.eachparticle(system)) - # Result - dv = zeros(3, 8) - TrixiParticles.interact!(dv, v, u, v, u, nhs, system, system) + # Result + dv = zeros(n_vars, n_particles) + TrixiParticles.interact!(dv, v, u, v, u, nhs, system, system) - # Linear momentum conservation - # ∑ m_a dv_a - deriv_linear_momentum = sum(mass[i]' .* view(dv, 1:2, :), dims=2) + # Linear momentum conservation + # ∑ m_a dv_a + deriv_linear_momentum = sum(fluid.mass' .* view(dv, 1:2, :), dims=2) - @test isapprox(deriv_linear_momentum, zeros(2, 1), atol=3e-14) + @test isapprox(deriv_linear_momentum, zeros(2, 1), atol=4e-14) - # Angular momentum conservation - # m_a (r_a × dv_a) - function deriv_angular_momentum(particle) - r_a = SVector(u[1, particle], u[2, particle], 0.0) - dv_a = SVector(dv[1, particle], dv[2, particle], 0.0) + # Angular momentum conservation + # m_a (r_a × dv_a) + function deriv_angular_momentum(particle) + r_a = SVector(u[1, particle], u[2, particle], 0.0) + dv_a = SVector(dv[1, particle], dv[2, particle], 0.0) - return mass[i][particle] * cross(r_a, dv_a) - end + return fluid.mass[particle] * cross(r_a, dv_a) + end - # ∑ m_a (r_a × dv_a) - deriv_angular_momentum = sum(deriv_angular_momentum, 1:8) + # ∑ m_a (r_a × dv_a) + deriv_angular_momentum = sum(deriv_angular_momentum, 1:n_particles) - @test isapprox(deriv_angular_momentum, zeros(3), atol=2e-14) + # Cross product is always 3-dimensional + @test isapprox(deriv_angular_momentum, zeros(3), atol=3e-15) - # Total energy conservation - drho(::ContinuityDensity, particle) = dv[3, particle] + # Total energy conservation + drho(::ContinuityDensity, particle) = dv[end, particle] - # Derivative of the density summation. This is a slightly different - # formulation of the continuity equation. - function drho_particle(particle, neighbor) - m_b = mass[i][neighbor] - vdiff = TrixiParticles.current_velocity(v, system, particle) - - TrixiParticles.current_velocity(v, system, neighbor) + function drho(::SummationDensity, particle) + return sum(neighbor -> drho_particle(particle, neighbor), 1:n_particles) + end - pos_diff = TrixiParticles.current_coords(u, system, particle) - - TrixiParticles.current_coords(u, system, neighbor) - distance = norm(pos_diff) + # Derivative of the density summation. This is a slightly different + # formulation of the continuity equation. + function drho_particle(particle, neighbor) + m_b = TrixiParticles.hydrodynamic_mass(system, neighbor) + v_diff = TrixiParticles.current_velocity(v, system, particle) - + TrixiParticles.current_velocity(v, system, neighbor) - # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return 0.0 + pos_diff = TrixiParticles.current_coords(u, system, particle) - + TrixiParticles.current_coords(u, system, neighbor) + distance = norm(pos_diff) - grad_kernel = TrixiParticles.smoothing_kernel_grad(system, - pos_diff, - distance) + # Only consider particles with a distance > 0 + distance < sqrt(eps()) && return 0.0 - return m_b * dot(vdiff, grad_kernel) - end + grad_kernel = TrixiParticles.smoothing_kernel_grad(system, pos_diff, + distance) - function drho(::SummationDensity, particle) - return sum(neighbor -> drho_particle(particle, neighbor), 1:8) - end + return m_b * dot(v_diff, grad_kernel) + end - # m_a (v_a ⋅ dv_a + dte_a), - # where `te` is the thermal energy, called `u` in the Price paper. - function deriv_energy(particle) - dte_a = pressure[i][particle] / density[i][particle]^2 * - drho(density_calculator, particle) - v_a = TrixiParticles.extract_svector(v, system, particle) - dv_a = TrixiParticles.extract_svector(dv, system, particle) + # m_a (v_a ⋅ dv_a + dte_a), + # where `te` is the thermal energy, called `u` in the Price paper. + function deriv_energy(particle) + p_a = fluid.pressure[particle] + rho_a = fluid.density[particle] + dte_a = p_a / rho_a^2 * drho(density_calculator, particle) + v_a = TrixiParticles.extract_svector(v, system, particle) + dv_a = TrixiParticles.extract_svector(dv, system, particle) - return mass[i][particle] * (dot(v_a, dv_a) + dte_a) - end + return fluid.mass[particle] * (dot(v_a, dv_a) + dte_a) + end - # ∑ m_a (v_a ⋅ dv_a + dte_a) - deriv_total_energy = sum(deriv_energy, 1:8) + # ∑ m_a (v_a ⋅ dv_a + dte_a) + deriv_total_energy = sum(deriv_energy, 1:n_particles) - @test isapprox(deriv_total_energy, 0.0, atol=4e-14) - end + @test isapprox(deriv_total_energy, 0.0, atol=2e-15) end end end diff --git a/test/setups/rectangular_tank.jl b/test/setups/rectangular_tank.jl index 278cc6951..8b9932101 100644 --- a/test/setups/rectangular_tank.jl +++ b/test/setups/rectangular_tank.jl @@ -10,6 +10,7 @@ @testset "Coordinates" begin particle_spacings = [0.1, 0.2] spacing_ratios = [1, 3] + min_coordinates = [(0.0, 0.0), (-0.3, 2.0)] expected_fluid_coords = [ [0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45 0.05 0.15 0.25 0.35 0.45; @@ -33,17 +34,22 @@ ], ] - @testset "Particle Spacing: $(particle_spacings[i])" for i in eachindex(particle_spacings) - @testset "Spacing Ratio: $(spacing_ratios[j])" for j in eachindex(spacing_ratios) - tank = RectangularTank(particle_spacings[i], - (water_width, water_height), - (tank_width, tank_height), - water_density, - spacing_ratio=spacing_ratios[j]) - - @test isapprox(tank.fluid.coordinates, expected_fluid_coords[i]) - @test isapprox(tank.boundary.coordinates, - expected_bound_coords[i][j]) + @testset "Move Tank: $(min_coordinates[i])" for i in eachindex(min_coordinates) + @testset "Particle Spacing: $(particle_spacings[j])" for j in eachindex(particle_spacings) + @testset "Spacing Ratio: $(spacing_ratios[k])" for k in eachindex(spacing_ratios) + tank = RectangularTank(particle_spacings[j], + (water_width, water_height), + (tank_width, tank_height), + water_density, + spacing_ratio=spacing_ratios[k], + min_coordinates=min_coordinates[i]) + expected_fluid_coords_ = copy(expected_fluid_coords[j]) + expected_bound_coords_ = copy(expected_bound_coords[j][k]) + expected_fluid_coords_ .+= min_coordinates[i] + expected_bound_coords_ .+= min_coordinates[i] + @test isapprox(tank.fluid.coordinates, expected_fluid_coords_) + @test isapprox(tank.boundary.coordinates, expected_bound_coords_) + end end end end @@ -267,6 +273,7 @@ end @testset "Coordinates" begin particle_spacings = [0.1, 0.2] spacing_ratios = [1, 3] + min_coordinates = [(0.0, 0.0, 0.0), (-0.3, 2.0, -0.5)] expected_fluid_coords = [ [0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35 0.05 0.15 0.25 0.35; @@ -296,17 +303,23 @@ end ], ] - @testset "Particle Spacing: $(particle_spacings[i])" for i in eachindex(particle_spacings) - @testset "Spacing Ratio: $(spacing_ratios[j])" for j in eachindex(spacing_ratios) - tank = RectangularTank(particle_spacings[i], - (water_width, water_height, water_depth), - (tank_width, tank_height, tank_depth), - water_density, - spacing_ratio=spacing_ratios[j]) - - @test isapprox(tank.fluid.coordinates, expected_fluid_coords[i]) - @test isapprox(tank.boundary.coordinates, - expected_bound_coords[i][j]) + @testset "Move Tank: $(min_coordinates[i])" for i in eachindex(min_coordinates) + @testset "Particle Spacing: $(particle_spacings[j])" for j in eachindex(particle_spacings) + @testset "Spacing Ratio: $(spacing_ratios[k])" for k in eachindex(spacing_ratios) + tank = RectangularTank(particle_spacings[j], + (water_width, water_height, water_depth), + (tank_width, tank_height, tank_depth), + water_density, + spacing_ratio=spacing_ratios[k], + min_coordinates=min_coordinates[i]) + expected_fluid_coords_ = copy(expected_fluid_coords[j]) + expected_bound_coords_ = copy(expected_bound_coords[j][k]) + expected_fluid_coords_ .+= min_coordinates[i] + expected_bound_coords_ .+= min_coordinates[i] + + @test isapprox(tank.fluid.coordinates, expected_fluid_coords_) + @test isapprox(tank.boundary.coordinates, expected_bound_coords_) + end end end end diff --git a/test/setups/sphere_shape.jl b/test/setups/sphere_shape.jl index e39cb2151..7577ed1ab 100644 --- a/test/setups/sphere_shape.jl +++ b/test/setups/sphere_shape.jl @@ -16,6 +16,7 @@ "Hollow RoundSphere with One Layer", "Hollow RoundSphere with Multiple Layers", "Hollow RoundSphere Outwards", + "Hollow RoundSphere Segment", ] shapes = [ @@ -39,6 +40,9 @@ SphereShape(particle_spacing, radius, (-3.0, 4.0), 1000.0, sphere_type=RoundSphere(), n_layers=3, layer_outwards=true), + SphereShape(particle_spacing, radius, (3.0, -4), 1000.0, + sphere_type=RoundSphere(start_angle=pi, end_angle=1.5pi), + n_layers=3, layer_outwards=true), ] expected_coords = [ @@ -66,6 +70,8 @@ -4.0 -3.22865487 -2.8182307 -2.96076952 -3.58957583 -4.41042417 -5.03923048 -5.1817693 -4.77134513 -4.0 -3.23463314 -2.58578644 -2.15224093 -2.0 -2.15224093 -2.58578644 -3.23463314 -4.0 -4.76536686 -5.41421356 -5.84775907 -6.0 -5.84775907 -5.41421356 -4.76536686], [-0.2 -0.31341967 -0.64449011 -1.16638994 -1.83683796 -2.60151845 -3.39848155 -4.16316204 -4.83361006 -5.35550989 -5.68658033 -5.8 -5.68658033 -5.35550989 -4.83361006 -4.16316204 -3.39848155 -2.60151845 -1.83683796 -1.16638994 -0.64449011 -0.31341967 0.6 0.50974048 0.24348792 -0.18540666 -0.75543671 -1.43801854 -2.19892464 -3.0 -3.80107536 -4.56198146 -5.24456329 -5.81459334 -6.24348792 -6.50974048 -6.6 -6.50974048 -6.24348792 -5.81459334 -5.24456329 -4.56198146 -3.80107536 -3.0 -2.19892464 -1.43801854 -0.75543671 -0.18540666 0.24348792 0.50974048 1.4 1.32929019 1.11943343 0.77717469 0.31351445 -0.25664487 -0.91497789 -1.64032522 -2.40937363 -3.19740525 -3.97909211 -4.72931014 -5.42394672 -6.04067566 -6.55967478 -6.96426302 -7.24143659 -7.38228689 -7.38228689 -7.24143659 -6.96426302 -6.55967478 -6.04067566 -5.42394672 -4.72931014 -3.97909211 -3.19740525 -2.40937363 -1.64032522 -0.91497789 -0.25664487 0.31351445 0.77717469 1.11943343 1.32929019; 4.0 4.78885116 5.51379429 6.11609881 6.54696959 6.77150004 6.77150004 6.54696959 6.11609881 5.51379429 4.78885116 4.0 3.21114884 2.48620571 1.88390119 1.45303041 1.22849996 1.22849996 1.45303041 1.88390119 2.48620571 3.21114884 4.0 4.80107536 5.56198146 6.24456329 6.81459334 7.24348792 7.50974048 7.6 7.50974048 7.24348792 6.81459334 6.24456329 5.56198146 4.80107536 4.0 3.19892464 2.43801854 1.75543671 1.18540666 0.75651208 0.49025952 0.4 0.49025952 0.75651208 1.18540666 1.75543671 2.43801854 3.19892464 4.0 4.78565034 5.54604923 6.25675682 6.89493039 7.44005852 7.87462034 8.18464867 8.36017895 8.39556949 8.28968281 8.0459222 7.67212232 7.1802974 6.58625511 5.90908845 5.17056212 4.39441296 3.60558704 2.82943788 2.09091155 1.41374489 0.8197026 0.32787768 -0.0459222 -0.28968281 -0.39556949 -0.36017895 -0.18464867 0.12537966 0.55994148 1.10506961 1.74324318 2.45395077 3.21434966], + [0.20000000000000018 0.33704175437357 0.7347524157501475 1.354201293581075 2.134752415750147 2.9999999999999996 -0.5999999999999996 -0.5097404838545647 -0.2434879244487087 0.18540666311509257 0.7554367133085589 1.4380185391767903 2.1989246377572673 2.9999999999999996 -1.4000000000000004 -1.3331541132537161 -1.1346475314579978 -0.8105117766515302 -0.37059554972350384 0.17173451737922596 0.8000000000000016 1.4951113693670546 2.2359480182655065 2.999999999999999; + -3.9999999999999996 -4.865247584249852 -5.645798706418924 -6.2652475842498525 -6.662958245626429 -6.8 -3.9999999999999996 -4.801075362242731 -5.561981460823208 -6.24456328669144 -6.814593336884906 -7.243487924448708 -7.509740483854564 -7.6 -3.9999999999999996 -4.764051981734492 -5.504888630632941 -6.200000000000001 -6.828265482620774 -7.370595549723503 -7.810511776651531 -8.134647531457997 -8.333154113253716 -8.4], ] @testset "$(shape_names[i])" for i in eachindex(shape_names) 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 diff --git a/test/test_util.jl b/test/test_util.jl index 6897bad29..c087a4dea 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -3,6 +3,7 @@ using TrixiParticles using LinearAlgebra using Printf using QuadGK: quadgk +using Random: Random """ @trixi_testset "name of the testset" #= code to test #= @@ -32,3 +33,35 @@ macro trixi_testset(name, expr) nothing end end + +function perturb!(data, amplitude) + for i in eachindex(data) + # Perturbation in the interval (-amplitude, amplitude) + data[i] += 2amplitude * rand() - amplitude + end + + return data +end + +# Rectangular patch of particles, optionally with a perturbation in position and quantities +function rectangular_patch(particle_spacing, size; density=1000.0, pressure=0.0, seed=1, + perturbation_factor=1.0) + # Fixed seed to ensure reproducibility + Random.seed!(seed) + + # Center particle at the origin (assuming odd size) + min_corner = -particle_spacing / 2 .* size + ic = RectangularShape(particle_spacing, size, min_corner, density, pressure=pressure) + + perturb!(ic.mass, perturbation_factor * 0.1 * ic.mass[1]) + perturb!(ic.density, perturbation_factor * 0.1density) + perturb!(ic.pressure, perturbation_factor * 2000) + perturb!(ic.velocity, perturbation_factor * 0.5particle_spacing) + perturb!(ic.coordinates, perturbation_factor * 0.5particle_spacing) + + # Don't perturb center particle position + center_particle = ceil(Int, prod(size) / 2) + ic.coordinates[:, center_particle] .= 0.0 + + return ic +end