diff --git a/src/general/corrections.jl b/src/general/corrections.jl deleted file mode 100644 index dd66d4efd..000000000 --- a/src/general/corrections.jl +++ /dev/null @@ -1,439 +0,0 @@ -# Sorted in order of computational cost -@doc raw""" - AkinciFreeSurfaceCorrection(rho0) - -Free surface correction according to [Akinci et al. (2013)](@cite Akinci2013). -At a free surface, the mean density is typically lower than the reference density, -resulting in reduced surface tension and viscosity forces. -The free surface correction adjusts the viscosity, pressure, and surface tension forces -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. -""" -struct AkinciFreeSurfaceCorrection{ELTYPE} - rho0::ELTYPE - - function AkinciFreeSurfaceCorrection(rho0) - ELTYPE = eltype(rho0) - return new{ELTYPE}(rho0) - end -end - -# `rho_mean` is the mean density of the fluid, which is used to determine correction values near the free surface. -# Return a tuple `(viscosity_correction, pressure_correction, surface_tension_correction)` representing the correction terms. -@inline function free_surface_correction(correction::AkinciFreeSurfaceCorrection, - particle_system, rho_mean) - # Equation 4 in ref - k = correction.rho0 / rho_mean - - # Viscosity, pressure, surface_tension - return k, 1, k -end - -@inline function free_surface_correction(correction, particle_system, rho_mean) - return 1, 1, 1 -end - -@doc raw""" - ShepardKernelCorrection() - -Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation -to obtain a 0-th order accurate result, which was first proposed by [Li et al. (1996)](@cite Li1996). - -The kernel correction coefficient is determined by -```math -c(x) = \sum_{b=1} V_b W_b(x), -``` -where ``V_b = m_b / \rho_b`` is the volume of particle ``b``. - -This correction is applied with [`SummationDensity`](@ref) to correct the density and leads -to an improvement, especially at free surfaces. - -!!! note - - It is also referred to as "0th order correction". - - In 2D, we can expect an increase of about 5--6% in computation time. -""" -struct ShepardKernelCorrection end - -@doc raw""" - KernelCorrection() - -Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation -to obtain a 0-th order accurate result, which was first proposed by Li et al. -This can be further extended to obtain a kernel corrected gradient as shown by [Basa et al. (2008)](@cite Basa2008). - -The kernel correction coefficient is determined by -```math -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) - 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)}. -``` - -This correction can be applied with [`SummationDensity`](@ref) and -[`ContinuityDensity`](@ref), which leads to an improvement, especially at free surfaces. - -!!! note - - This only works when the boundary model uses [`SummationDensity`](@ref) (yet). - - It is also referred to as "0th order correction". - - In 2D, we can expect an increase of about 10--15% in computation time. -""" -struct KernelCorrection end - -@doc raw""" - MixedKernelGradientCorrection() - -Combines [`GradientCorrection`](@ref) and [`KernelCorrection`](@ref), -which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1999)). - -# Notes: -- Stability issues, especially when particles separate into small clusters. -- Doubles the computational effort. -""" -struct MixedKernelGradientCorrection end - -function kernel_correction_coefficient(system::AbstractFluidSystem, particle) - return system.cache.kernel_correction_coefficient[particle] -end - -function kernel_correction_coefficient(system::AbstractBoundarySystem, particle) - return system.boundary_model.cache.kernel_correction_coefficient[particle] -end - -function compute_correction_values!(system, correction, u, v_ode, u_ode, semi) - return system -end - -function compute_correction_values!(system, ::ShepardKernelCorrection, u, v_ode, u_ode, - semi) - return compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, - semi, - system.cache.kernel_correction_coefficient) -end - -function compute_correction_values!(system::AbstractBoundarySystem, - ::ShepardKernelCorrection, u, - v_ode, u_ode, semi) - return compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, - semi, - system.boundary_model.cache.kernel_correction_coefficient) -end - -function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi, - kernel_correction_coefficient) - set_zero!(kernel_correction_coefficient) - - # Use all other systems for the density summation - @trixi_timeit timer() "compute correction value" 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) - - # Loop over all pairs of particles and neighbors within the kernel cutoff - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi) do particle, neighbor, pos_diff, distance - rho_b = current_density(v_neighbor_system, neighbor_system, neighbor) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - volume = m_b / rho_b - - kernel_correction_coefficient[particle] += volume * - smoothing_kernel(system, distance, - particle) - end - end - - return kernel_correction_coefficient -end - -function dw_gamma(system::AbstractFluidSystem, particle) - return extract_svector(system.cache.dw_gamma, system, particle) -end - -function dw_gamma(system::AbstractBoundarySystem, particle) - return extract_svector(system.boundary_model.cache.dw_gamma, system, particle) -end - -function compute_correction_values!(system::AbstractFluidSystem, - correction::Union{KernelCorrection, - MixedKernelGradientCorrection}, u, - v_ode, u_ode, semi) - compute_correction_values!(system, correction, current_coordinates(u, system), v_ode, - u_ode, semi, - system.cache.kernel_correction_coefficient, - system.cache.dw_gamma) -end - -function compute_correction_values!(system::AbstractBoundarySystem, - correction::Union{KernelCorrection, - MixedKernelGradientCorrection}, u, - v_ode, u_ode, semi) - compute_correction_values!(system, correction, current_coordinates(u, system), v_ode, - u_ode, semi, - system.boundary_model.cache.kernel_correction_coefficient, - system.boundary_model.cache.dw_gamma) -end - -function compute_correction_values!(system, - ::Union{KernelCorrection, - MixedKernelGradientCorrection}, system_coords, - v_ode, - u_ode, semi, kernel_correction_coefficient, dw_gamma) - set_zero!(kernel_correction_coefficient) - set_zero!(dw_gamma) - - # Use all other systems for the density summation - @trixi_timeit timer() "compute correction value" 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) - - # Loop over all pairs of particles and neighbors within the kernel cutoff - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi) do particle, neighbor, pos_diff, distance - rho_b = current_density(v_neighbor_system, neighbor_system, neighbor) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - volume = m_b / rho_b - - # Use uncorrected kernel to compute correction coefficients - W = kernel(system_smoothing_kernel(system), distance, - smoothing_length(system, particle)) - - kernel_correction_coefficient[particle] += volume * W - - # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details. - if distance^2 > eps(initial_smoothing_length(system)^2) - grad_W = kernel_grad(system_smoothing_kernel(system), pos_diff, distance, - smoothing_length(system, particle)) - tmp = volume * grad_W - for i in axes(dw_gamma, 1) - dw_gamma[i, particle] += tmp[i] - end - end - end - end - - for particle in eachparticle(system), i in axes(dw_gamma, 1) - dw_gamma[i, particle] /= kernel_correction_coefficient[particle] - end -end - -@doc raw""" - GradientCorrection() - -Compute the corrected gradient of particle interactions based on their relative positions -(see [Bonet, 1999](@cite Bonet1999)). - -# Mathematical Details - -Given the standard SPH representation, the gradient of a field ``A`` at particle ``a`` is -given by - -```math -\nabla A_a = \sum_b m_b \frac{A_b - A_a}{\rho_b} \nabla_{r_a} W(\Vert r_a - r_b \Vert, h), -``` -where ``m_b`` is the mass of particle ``b`` and ``\rho_b`` is the density of particle ``b``. - -The gradient correction, as commonly proposed, involves multiplying this gradient with a correction matrix $L$: - -```math -\tilde{\nabla} A_a = \bm{L}_a \nabla A_a -``` - -The correction matrix $\bm{L}_a$ is computed based on the provided particle configuration, -aiming to make the corrected gradient more accurate, especially near domain boundaries. - -To satisfy -```math -\sum_b V_b r_{ba} \otimes \tilde{\nabla}W_b(r_a) = \left( \sum_b V_b r_{ba} \otimes \nabla W_b(r_a) \right) \bm{L}_a^T = \bm{I} -``` -the correction matrix $\bm{L}_a$ is evaluated explicitly as -```math -\bm{L}_a = \left( \sum_b V_b \nabla W_b(r_{a}) \otimes r_{ba} \right)^{-1}. -``` - -!!! 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. -""" -struct GradientCorrection end - -@doc raw""" - BlendedGradientCorrection() - -Calculate a blended gradient to reduce the stability issues of the [`GradientCorrection`](@ref) -as explained by [Bonet (1999)](@cite Bonet1999). - -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, system, coordinates, density_fun, - semi) - (; mass) = system - - set_zero!(corr_matrix) - - # Loop over all pairs of particles and neighbors within the kernel cutoff. - foreach_point_neighbor(system, system, coordinates, coordinates, - semi) do particle, neighbor, pos_diff, distance - volume = mass[neighbor] / density_fun(neighbor) - - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - iszero(grad_kernel) && return - - result = volume * grad_kernel * pos_diff' - - @inbounds for j in 1:ndims(system), i in 1:ndims(system) - corr_matrix[i, j, particle] -= result[i, j] - end - end - - correction_matrix_inversion_step!(corr_matrix, system, semi) - - return corr_matrix -end - -function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, - coordinates, v_ode, u_ode, semi, - correction, 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) - - foreach_point_neighbor(system, neighbor_system, coordinates, neighbor_coords, - semi) do particle, neighbor, pos_diff, distance - volume = hydrodynamic_mass(neighbor_system, neighbor) / - current_density(v_neighbor_system, neighbor_system, neighbor) - smoothing_length_ = smoothing_length(system, particle) - - function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, - smoothing_length_, system, particle) - return smoothing_kernel_grad(system, pos_diff, distance, particle) - 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, semi) - - return corr_matrix -end - -function correction_matrix_inversion_step!(corr_matrix, system, semi) - @threaded semi for particle in eachparticle(system) - L = extract_smatrix(corr_matrix, system, particle) - - # The matrix `L` only becomes singular when the particle and all neighbors - # are collinear (in 2D) or lie all in the same plane (in 3D). - # This happens only when two (in 2D) or three (in 3D) particles are isolated, - # or in cases where there is only one layer of fluid particles on a wall. - # In these edge cases, we just disable the correction and set the corrected - # gradient to be the uncorrected one by setting `L` to the identity matrix. - # - # Proof: `L` is just a sum of tensor products of relative positions X_ab with - # themselves. According to - # https://en.wikipedia.org/wiki/Outer_product#Connection_with_the_matrix_product - # the sum of tensor products can be rewritten as A A^T, where the columns of A - # are the relative positions X_ab. The rank of A A^T is equal to the rank of A, - # so `L` is singular if and only if the position vectors X_ab don't span the - # full space, i.e., particle a and all neighbors lie on the same line (in 2D) - # or plane (in 3D). - if abs(det(L)) < 1.0f-9 - L_inv = I - else - L_inv = inv(L) - end - - # Write inverse back to `corr_matrix` - for j in 1:ndims(system), i in 1:ndims(system) - @inbounds corr_matrix[i, j, particle] = L_inv[i, j] - end - end - - return corr_matrix -end - -create_cache_correction(correction, density, NDIMS, nparticles) = (;) - -function create_cache_correction(::ShepardKernelCorrection, density, NDIMS, n_particles) - return (; kernel_correction_coefficient=similar(density)) -end - -function create_cache_correction(::KernelCorrection, density, NDIMS, n_particles) - dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles) - return (; kernel_correction_coefficient=similar(density), dw_gamma) -end - -function create_cache_correction(::Union{GradientCorrection, BlendedGradientCorrection}, - density, - NDIMS, n_particles) - correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles) - return (; correction_matrix) -end - -function create_cache_correction(::MixedKernelGradientCorrection, density, NDIMS, - n_particles) - dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles) - correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles) - - return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) -end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl deleted file mode 100644 index ceba8b5bd..000000000 --- a/src/general/semidiscretization.jl +++ /dev/null @@ -1,1121 +0,0 @@ -""" - Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch{NDIMS}()) - -The semidiscretization couples the passed systems to one simulation. - -# Arguments -- `systems`: Systems to be coupled in this semidiscretization - -# Keywords -- `neighborhood_search`: The neighborhood search to be used in the simulation. - By default, the [`GridNeighborhoodSearch`](@ref) is used. - Use `nothing` to loop over all particles (no neighborhood search). - To use other neighborhood search implementations, pass a template - of a neighborhood search. See [`copy_neighborhood_search`](@ref) - and the examples below for more details. - To use a periodic domain, pass a [`PeriodicBox`](@ref) to the - neighborhood search. -- `threaded_nhs_update=true`: Can be used to deactivate thread parallelization in the neighborhood search update. - This can be one of the largest sources of variations between simulations - with different thread numbers due to particle ordering changes. - -# Examples -```jldoctest; output = false, setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); ref_system = fluid_system) -semi = Semidiscretization(fluid_system, boundary_system) - -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=SerialUpdate())) - -periodic_box = PeriodicBox(min_corner = [0.0, 0.0], max_corner = [1.0, 1.0]) -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box)) - -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=PrecomputedNeighborhoodSearch{2}()) - -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=nothing) - -# output -┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ -│ Semidiscretization │ -│ ══════════════════ │ -│ #spatial dimensions: ………………………… 2 │ -│ #systems: ……………………………………………………… 2 │ -│ neighborhood search: ………………………… TrivialNeighborhoodSearch │ -│ total #particles: ………………………………… 636 │ -└──────────────────────────────────────────────────────────────────────────────────────────────────┘ -``` -""" -struct Semidiscretization{BACKEND, S, RU, RV, NS, UCU, IT} - systems :: S - ranges_u :: RU - ranges_v :: RV - neighborhood_searches :: NS - parallelization_backend :: BACKEND - update_callback_used :: UCU - integrate_tlsph :: IT # `false` if TLSPH integration is decoupled - - # Dispatch at `systems` to distinguish this constructor from the one below when - # 4 systems are passed. - # This is an internal constructor only used in `test/count_allocations.jl` - # and by Adapt.jl. - function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches, - parallelization_backend::PointNeighbors.ParallelizationBackend, - update_callback_used, integrate_tlsph) - new{typeof(parallelization_backend), typeof(systems), typeof(ranges_u), - typeof(ranges_v), typeof(neighborhood_searches), - typeof(update_callback_used), - typeof(integrate_tlsph)}(systems, ranges_u, ranges_v, - neighborhood_searches, parallelization_backend, - update_callback_used, integrate_tlsph) - end -end - -function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; - neighborhood_search=GridNeighborhoodSearch{ndims(first(systems))}(), - parallelization_backend=PolyesterBackend()) - systems = filter(system -> !isnothing(system), systems) - - # Check e.g. that the boundary systems are using a state equation if EDAC is not used. - # Other checks might be added here later. - check_configuration(systems, neighborhood_search) - - sizes_u = [u_nvariables(system) * n_integrated_particles(system) - for system in systems] - ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) - for i in eachindex(sizes_u)) - sizes_v = [v_nvariables(system) * n_integrated_particles(system) - for system in systems] - ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) - for i in eachindex(sizes_v)) - - # Create a tuple of n neighborhood searches for each of the n systems. - # We will need one neighborhood search for each pair of systems. - searches = Tuple(Tuple(create_neighborhood_search(neighborhood_search, - system, neighbor) - for neighbor in systems) - for system in systems) - - # These will be set to true inside the `UpdateCallback`. - # Some techniques require the use of this callback, and this flag can be used - # to determine if the callback is used in a simulation. - update_callback_used = Ref(false) - - # Always integrate TLSPH systems together with other systems. - # For split integration, a copy of the semidiscretization will be created - # with this set to false. - integrate_tlsph = Ref(true) - - return Semidiscretization(systems, ranges_u, ranges_v, searches, - parallelization_backend, update_callback_used, - integrate_tlsph) -end - -# Inline show function e.g. Semidiscretization(neighborhood_search=...) -function Base.show(io::IO, semi::Semidiscretization) - @nospecialize semi # reduce precompilation time - - print(io, "Semidiscretization(") - for system in semi.systems - print(io, system, ", ") - end - print(io, "neighborhood_search=") - print(io, semi.neighborhood_searches |> eltype |> eltype |> nameof) - print(io, ")") -end - -# Show used during summary printout -function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) - @nospecialize semi # reduce precompilation time - - if get(io, :compact, false) - show(io, semi) - else - summary_header(io, "Semidiscretization") - summary_line(io, "#spatial dimensions", ndims(semi.systems[1])) - summary_line(io, "#systems", length(semi.systems)) - summary_line(io, "neighborhood search", - semi.neighborhood_searches |> eltype |> eltype |> nameof) - summary_line(io, "total #particles", sum(nparticles.(semi.systems))) - summary_footer(io) - end -end - -function create_neighborhood_search(::Nothing, system, neighbor) - nhs = TrivialNeighborhoodSearch{ndims(system)}() - - return create_neighborhood_search(nhs, system, neighbor) -end - -function create_neighborhood_search(neighborhood_search, system, neighbor) - return copy_neighborhood_search(neighborhood_search, compact_support(system, neighbor), - nparticles(neighbor)) -end - -@inline function compact_support(system, neighbor) - (; smoothing_kernel) = system - # TODO: Variable search radius for NHS? - return compact_support(smoothing_kernel, initial_smoothing_length(system)) -end - -@inline function compact_support(system::OpenBoundarySystem, - neighbor::OpenBoundarySystem) - # This NHS is never used - return zero(eltype(system)) -end - -@inline function compact_support(system::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}, - neighbor::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}) - # Use the compact support of the fluid - return compact_support(system.fluid_system, neighbor.fluid_system) -end - -@inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem) - # This NHS is never used - return zero(eltype(system)) -end - -@inline function compact_support(system::BoundaryDEMSystem, neighbor::DEMSystem) - # Use the compact support of the DEMSystem - return compact_support(neighbor, system) -end - -@inline function compact_support(system::TotalLagrangianSPHSystem, - neighbor::TotalLagrangianSPHSystem) - (; smoothing_kernel, smoothing_length) = system - return compact_support(smoothing_kernel, smoothing_length) -end - -@inline function compact_support(system::Union{TotalLagrangianSPHSystem, - WallBoundarySystem}, - neighbor) - return compact_support(system, system.boundary_model, neighbor) -end - -@inline function compact_support(system, model::BoundaryModelMonaghanKajtar, neighbor) - # Use the compact support of the fluid for structure-fluid interaction - return compact_support(neighbor, system) -end - -@inline function compact_support(system, model::BoundaryModelMonaghanKajtar, - neighbor::WallBoundarySystem) - # This NHS is never used - return zero(eltype(system)) -end - -@inline function compact_support(system, model::BoundaryModelDummyParticles, neighbor) - # TODO: Monaghan-Kajtar BC are using the fluid's compact support for structure-fluid - # interaction. Dummy particle BC use the model's compact support, which is also used - # for density summations. - (; smoothing_kernel, smoothing_length) = model - return compact_support(smoothing_kernel, smoothing_length) -end - -@inline function get_neighborhood_search(system, semi) - (; neighborhood_searches) = semi - - system_index = system_indices(system, semi) - - return neighborhood_searches[system_index][system_index] -end - -@inline function get_neighborhood_search(system, neighbor_system, semi) - (; neighborhood_searches) = semi - - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor_system, semi) - - return neighborhood_searches[system_index][neighbor_index] -end - -@inline function system_indices(system, semi) - # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict` - # is ~30x slower because `hash(::System)` is very slow. - index = findfirst(==(system), semi.systems) - - if isnothing(index) - throw(ArgumentError("system is not in the semidiscretization")) - end - - return index -end - -# This is just for readability to loop over all systems without allocations -@inline function foreach_system(f, semi::Union{NamedTuple, Semidiscretization}) - return foreach_noalloc(f, semi.systems) -end - -@inline foreach_system(f, systems) = foreach_noalloc(f, systems) - -""" - semidiscretize(semi, tspan; reset_threads=true) - -Create an `ODEProblem` from the semidiscretization with the specified `tspan`. - -# Arguments -- `semi`: A [`Semidiscretization`](@ref) holding the systems involved in the simulation. -- `tspan`: The time span over which the simulation will be run. - -# Keywords -- `reset_threads`: A boolean flag to reset Polyester.jl threads before the simulation (default: `true`). - After an error within a threaded loop, threading might be disabled. Resetting the threads before the simulation - ensures that threading is enabled again for the simulation. - See also [trixi-framework/Trixi.jl#1583](https://github.com/trixi-framework/Trixi.jl/issues/1583). - -# Returns -A `DynamicalODEProblem` (see [the OrdinaryDiffEq.jl docs](https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/)) -to be integrated with [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl). -Note that this is not a true `DynamicalODEProblem` where the acceleration does not depend on the velocity. -Therefore, not all integrators designed for `DynamicalODEProblem`s will work properly. -However, all integrators designed for `ODEProblem`s can be used. -See [time integration](@ref time_integration) for more details. - -# Examples -```jldoctest; output = false, filter = r"u0: .*", setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); ref_system = fluid_system) -semi = Semidiscretization(fluid_system, boundary_system) -tspan = (0.0, 1.0) -ode_problem = semidiscretize(semi, tspan) - -# output -ODEProblem with uType RecursiveArrayTools.ArrayPartition{Float64, Tuple{TrixiParticles.ThreadedBroadcastArray{Float64, 1, Vector{Float64}, PolyesterBackend}, TrixiParticles.ThreadedBroadcastArray{Float64, 1, Vector{Float64}, PolyesterBackend}}} and tType Float64. In-place: true -Non-trivial mass matrix: false -timespan: (0.0, 1.0) -u0: ([...], [...]) *this line is ignored by filter* -``` -""" -function semidiscretize(semi, tspan; reset_threads=true) - (; systems) = semi - - @assert all(system -> eltype(system) === eltype(systems[1]), systems) - ELTYPE = eltype(systems[1]) - - # Optionally reset Polyester.jl threads. See - # https://github.com/trixi-framework/Trixi.jl/issues/1583 - # https://github.com/JuliaSIMD/Polyester.jl/issues/30 - if reset_threads - Polyester.reset_threads!() - end - - sizes_u = (u_nvariables(system) * n_integrated_particles(system) for system in systems) - sizes_v = (v_nvariables(system) * n_integrated_particles(system) for system in systems) - - # Use either the specified backend, e.g., `CUDABackend` or `MetalBackend` or - # use CPU vectors for all CPU backends. - u0_ode_ = allocate(semi.parallelization_backend, ELTYPE, sum(sizes_u)) - v0_ode_ = allocate(semi.parallelization_backend, ELTYPE, sum(sizes_v)) - - if semi.parallelization_backend isa KernelAbstractions.Backend - u0_ode = u0_ode_ - v0_ode = v0_ode_ - else - # CPU vectors are wrapped in `ThreadedBroadcastArray`s - # to make broadcasting (which is done by OrdinaryDiffEq.jl) multithreaded. - # See https://github.com/trixi-framework/TrixiParticles.jl/pull/722 for more details. - u0_ode = ThreadedBroadcastArray(u0_ode_; - parallelization_backend=semi.parallelization_backend) - v0_ode = ThreadedBroadcastArray(v0_ode_; - parallelization_backend=semi.parallelization_backend) - end - - # Set initial condition - foreach_system(semi) do system - u0_system = wrap_u(u0_ode, system, semi) - v0_system = wrap_v(v0_ode, system, semi) - - write_u0!(u0_system, system) - write_v0!(v0_system, system) - end - - # TODO initialize after adapting to the GPU. - # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86. - initialize_neighborhood_searches!(semi) - - if semi.parallelization_backend isa KernelAbstractions.Backend - # Convert all arrays to the correct array type. - # When e.g. `parallelization_backend=CUDABackend()`, this will convert all `Array`s - # to `CuArray`s, moving data to the GPU. - # See the comments in general/gpu.jl for more details. - semi_ = Adapt.adapt(semi.parallelization_backend, semi) - - # We now have a new `Semidiscretization` with new systems. - # This means that systems linking to other systems still point to old systems. - # Therefore, we have to re-link them, which yields yet another `Semidiscretization`. - # Note that this re-creates systems containing links, so it only works as long - # as systems don't link to other systems containing links. - semi_new = Semidiscretization(set_system_links.(semi_.systems, Ref(semi_)), - semi_.ranges_u, semi_.ranges_v, - semi_.neighborhood_searches, - semi_.parallelization_backend, - semi_.update_callback_used, semi_.integrate_tlsph) - else - semi_new = semi - end - - # Initialize all particle systems - foreach_system(semi_new) do system - # Initialize this system - initialize!(system, semi_new) - end - - # Reset callback flag that will be set by the `UpdateCallback` - semi_new.update_callback_used[] = false - - return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new) -end - -""" - restart_with!(semi, sol) - -Set the initial coordinates and velocities of all systems in `semi` to the final values -in the solution `sol`. -[`semidiscretize`](@ref) has to be called again afterwards, or another -[`Semidiscretization`](@ref) can be created with the updated systems. - -# Arguments -- `semi`: The semidiscretization -- `sol`: The `ODESolution` returned by `solve` of `OrdinaryDiffEq` -""" -function restart_with!(semi, sol; reset_threads=true) - # Optionally reset Polyester.jl threads. See - # https://github.com/trixi-framework/Trixi.jl/issues/1583 - # https://github.com/JuliaSIMD/Polyester.jl/issues/30 - if reset_threads - Polyester.reset_threads!() - end - - initialize_neighborhood_searches!(semi) - - foreach_system(semi) do system - v = wrap_v(sol.u[end].x[1], system, semi) - u = wrap_u(sol.u[end].x[2], system, semi) - - restart_with!(system, v, u) - end - - # Reset callback flag that will be set by the `UpdateCallback` - semi.update_callback_used[] = false - - return semi -end - -function initialize_neighborhood_searches!(semi) - foreach_system(semi) do system - foreach_system(semi) do neighbor - # TODO Initialize after adapting to the GPU. - # Currently, this cannot use `semi.parallelization_backend` - # because data is still on the CPU. - PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi), - initial_coordinates(system), - initial_coordinates(neighbor), - eachindex_y=each_active_particle(neighbor), - parallelization_backend=PolyesterBackend()) - end - end - - return semi -end - -# We have to pass `system` here for type stability, -# since the type of `system` determines the return type. -@inline function wrap_v(v_ode, system, semi) - (; ranges_v) = semi - - range = ranges_v[system_indices(system, semi)] - - @boundscheck @assert length(range) == - v_nvariables(system) * n_integrated_particles(system) - - return wrap_array(v_ode, range, - (StaticInt(v_nvariables(system)), n_integrated_particles(system))) -end - -@inline function wrap_u(u_ode, system, semi) - (; ranges_u) = semi - - range = ranges_u[system_indices(system, semi)] - - @boundscheck @assert length(range) == - u_nvariables(system) * n_integrated_particles(system) - - return wrap_array(u_ode, range, - (StaticInt(u_nvariables(system)), n_integrated_particles(system))) -end - -@inline function wrap_array(array::Array, range, size) - # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(array), 2}, pointer(view(array, range)), size) - return PtrArray(pointer(view(array, range)), size) -end - -@inline function wrap_array(array::ThreadedBroadcastArray, range, size) - return ThreadedBroadcastArray(wrap_array(parent(array), range, size)) -end - -@inline function wrap_array(array, range, size) - # For non-`Array`s (typically GPU arrays), just reshape. Calling the `PtrArray` code - # above for a `CuArray` yields another `CuArray` (instead of a `PtrArray`) - # and is 8 times slower with double the allocations. - # - # Note that `size` might contain `StaticInt`s, so convert to `Int` first. - return reshape(view(array, range), Int.(size)) -end - -function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) - (; systems) = semi - - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system, semi), systems) -end - -function drift!(du_ode, v_ode, u_ode, semi, t) - @trixi_timeit timer() "drift!" begin - @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode) - - @trixi_timeit timer() "velocity" begin - # Set velocity and add acceleration for each system - foreach_system(semi) do system - du = wrap_u(du_ode, system, semi) - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - @threaded semi for particle in each_integrated_particle(system) - # This can be dispatched per system - add_velocity!(du, v, u, particle, system, semi, t) - end - end - end - end - - return du_ode -end - -@inline function add_velocity!(du, v, u, particle, system, semi::Semidiscretization, t) - add_velocity!(du, v, u, particle, system, t) -end - -@inline function add_velocity!(du, v, u, particle, system::TotalLagrangianSPHSystem, - semi::Semidiscretization, t) - # Only add velocity for TLSPH systems if they are integrated - if semi.integrate_tlsph[] - add_velocity!(du, v, u, particle, system, t) - end -end - -@inline function add_velocity!(du, v, u, particle, system, t) - # Generic fallback for all systems that don't define this function - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] - end - - return du -end - -# Solid wall boundary system doesn't integrate the particle positions -@inline add_velocity!(du, v, u, particle, system::WallBoundarySystem, t) = du - -@inline function add_velocity!(du, v, u, particle, system::AbstractFluidSystem, t) - # This is zero unless a shifting technique is used - delta_v_ = delta_v(system, particle) - - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] + delta_v_[i] - end - - return du -end - -function kick!(dv_ode, v_ode, u_ode, semi, t) - @trixi_timeit timer() "kick!" begin - # Check that the `UpdateCallback` is used if required - check_update_callback(semi) - - @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode) - - @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, - semi, t) - - @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, - semi) - - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, - semi, t) - end - - return dv_ode -end - -# Update the systems and neighborhood searches (NHS) for a simulation -# before calling `interact!` to compute forces. -function update_systems_and_nhs(v_ode, u_ode, semi, t) - # First update step before updating the NHS - # (for example for writing the current coordinates in the TLSPH system) - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_positions!(system, v, u, v_ode, u_ode, semi, t) - end - - # Update NHS - @trixi_timeit timer() "update nhs" update_nhs!(semi, u_ode) - - # Second update step. - # This is used to calculate density and pressure of the fluid systems - # before updating the boundary systems, - # since the fluid pressure is needed by the Adami interpolation. - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_quantities!(system, v, u, v_ode, u_ode, semi, t) - end - - update_implicit_sph!(semi, v_ode, u_ode, t) - - # Perform correction and pressure calculation - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_pressure!(system, v, u, v_ode, u_ode, semi, t) - end - - # This update depends on the computed quantities of the fluid system and therefore - # needs to be after `update_quantities!`. - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t) - end - - # Final update step for all remaining systems - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_final!(system, v, u, v_ode, u_ode, semi, t) - end -end - -function update_nhs!(semi, u_ode) - # Update NHS for each pair of systems - foreach_system(semi) do system - u_system = wrap_u(u_ode, system, semi) - - foreach_system(semi) do neighbor - u_neighbor = wrap_u(u_ode, neighbor, semi) - neighborhood_search = get_neighborhood_search(system, neighbor, semi) - - update_nhs!(neighborhood_search, system, neighbor, u_system, u_neighbor, semi) - end - end -end - -# The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different -# semidiscretization for wrapping arrays. -# TODO `semi` is not used yet, but will be used when the source terms API is modified -# to match the custom quantities API. -function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) - foreach_system(semi_wrap) do system - dv = wrap_v(dv_ode, system, semi_wrap) - v = wrap_v(v_ode, system, semi_wrap) - u = wrap_u(u_ode, system, semi_wrap) - - @threaded semi for particle in each_integrated_particle(system) - # Dispatch by system type to exclude boundary systems. - # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function - # can be used in the `SplitIntegrationCallback` as well. - add_acceleration!(dv, particle, system, semi_wrap.integrate_tlsph[]) - add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t, - semi_wrap.integrate_tlsph[]) - end - end - - return dv_ode -end - -@inline source_terms(system) = nothing -@inline source_terms(system::Union{AbstractFluidSystem, AbstractStructureSystem}) = system.source_terms - -@inline function add_acceleration!(dv, particle, system, integrate_tlsph) - add_acceleration!(dv, particle, system) -end - -@inline function add_acceleration!(dv, particle, system::TotalLagrangianSPHSystem, - integrate_tlsph) - integrate_tlsph && add_acceleration!(dv, particle, system) -end - -@inline add_acceleration!(dv, particle, system) = dv - -@inline function add_acceleration!(dv, particle, - system::Union{AbstractFluidSystem, - AbstractStructureSystem}) - (; acceleration) = system - - for i in 1:ndims(system) - dv[i, particle] += acceleration[i] - end - - return dv -end - -@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t, - integrate_tlsph) - add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) -end - -@inline function add_source_terms_inner!(dv, v, u, particle, - system::TotalLagrangianSPHSystem, - source_terms_, t, integrate_tlsph) - integrate_tlsph && add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) -end - -@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) - coords = current_coords(u, system, particle) - velocity = current_velocity(v, system, particle) - density = current_density(v, system, particle) - pressure = current_pressure(v, system, particle) - - source = source_terms_(coords, velocity, density, pressure, t) - - # Loop over `eachindex(source)`, so that users could also pass source terms for - # the density when using `ContinuityDensity`. - for i in eachindex(source) - dv[i, particle] += source[i] - end - - return dv -end - -@inline add_source_terms_inner!(dv, v, u, particle, system, source_terms_::Nothing, t) = dv - -@doc raw""" - SourceTermDamping(; damping_coefficient) - -A source term to be used when a damping step is required before running a full simulation. -The term ``-c \cdot v_a`` is added to the acceleration ``\frac{\mathrm{d}v_a}{\mathrm{d}t}`` -of particle ``a``, where ``c`` is the damping coefficient and ``v_a`` is the velocity of -particle ``a``. - -# Keywords -- `damping_coefficient`: The coefficient ``d`` above. A higher coefficient means more - damping. A coefficient of `1e-4` is a good starting point for - damping a fluid at rest. - -# Examples -```jldoctest; output = false -source_terms = SourceTermDamping(; damping_coefficient=1e-4) - -# output -SourceTermDamping{Float64}(0.0001) -``` -""" -struct SourceTermDamping{ELTYPE} - damping_coefficient::ELTYPE - - function SourceTermDamping(; damping_coefficient) - return new{typeof(damping_coefficient)}(damping_coefficient) - end -end - -@inline function (source_term::SourceTermDamping)(coords, velocity, density, pressure, t) - (; damping_coefficient) = source_term - - return -damping_coefficient * velocity -end - -function system_interaction!(dv_ode, v_ode, u_ode, semi) - # Call `interact!` for each pair of systems - foreach_system(semi) do system - foreach_system(semi) do neighbor - # Construct string for the interactions timer. - # Avoid allocations from string construction when no timers are used. - if timeit_debug_enabled() - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor, semi) - timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" - else - timer_str = "" - end - - interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) - end - end - - return dv_ode -end - -# Function barrier to make benchmarking interactions easier. -# One can benchmark, e.g. the fluid-fluid interaction, with: -# dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x; -# @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi); -@inline function interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; timer_str="") - dv = wrap_v(dv_ode, system, semi) - v_system = wrap_v(v_ode, system, semi) - u_system = wrap_u(u_ode, system, semi) - - v_neighbor = wrap_v(v_ode, neighbor, semi) - u_neighbor = wrap_u(u_ode, neighbor, semi) - - @trixi_timeit timer() timer_str begin - interact!(dv, v_system, u_system, v_neighbor, u_neighbor, system, neighbor, semi) - end -end - -# NHS updates -# To prevent hard-to-find bugs, there is not default version -function update_nhs!(neighborhood_search, - system::AbstractFluidSystem, - neighbor::Union{AbstractFluidSystem, TotalLagrangianSPHSystem}, - u_system, u_neighbor, semi) - # The current coordinates of fluids and structures change over time - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=each_active_particle(neighbor)) -end - -function update_nhs!(neighborhood_search, - system::Union{AbstractFluidSystem, - OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - neighbor::WallBoundarySystem, - u_system, u_neighbor, semi) - # Boundary coordinates only change over time when `neighbor.ismoving[]` - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, neighbor.ismoving[])) -end - -function update_nhs!(neighborhood_search, - system::AbstractFluidSystem, neighbor::OpenBoundarySystem, - u_system, u_neighbor, semi) - # The current coordinates of fluids and open boundaries change over time. - - # TODO: Update only `active_coordinates` of open boundaries. - # Problem: Removing inactive particles from neighboring lists is necessary. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=each_active_particle(neighbor)) -end - -function update_nhs!(neighborhood_search, - system::OpenBoundarySystem, neighbor::AbstractFluidSystem, - u_system, u_neighbor, semi) - # The current coordinates of both open boundaries and fluids change over time. - - # TODO: Update only `active_coordinates` of open boundaries. - # Problem: Removing inactive particles from neighboring lists is necessary. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=each_active_particle(neighbor)) -end - -function update_nhs!(neighborhood_search, - system::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}, - neighbor::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}, - u_system, u_neighbor, semi) - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=each_active_particle(neighbor)) -end - -function update_nhs!(neighborhood_search, - system::OpenBoundarySystem, neighbor::TotalLagrangianSPHSystem, - u_system, u_neighbor, semi) - # Don't update. This NHS is never used. - return neighborhood_search -end - -function update_nhs!(neighborhood_search, - system::TotalLagrangianSPHSystem, neighbor::OpenBoundarySystem, - u_system, u_neighbor, semi) - # Don't update. This NHS is never used. - return neighborhood_search -end - -function update_nhs!(neighborhood_search, - system::TotalLagrangianSPHSystem, neighbor::AbstractFluidSystem, - u_system, u_neighbor, semi) - # The current coordinates of fluids and structured change over time - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=each_active_particle(neighbor)) -end - -function update_nhs!(neighborhood_search, - system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem, - u_system, u_neighbor, semi) - # Don't update. Neighborhood search works on the initial coordinates, which don't change. - return neighborhood_search -end - -function update_nhs!(neighborhood_search, - system::TotalLagrangianSPHSystem, neighbor::WallBoundarySystem, - u_system, u_neighbor, semi) - # The current coordinates of structured change over time. - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, neighbor.ismoving[])) -end - -# This function is the same as the one below to avoid ambiguous dispatch when using `Union` -function update_nhs!(neighborhood_search, - system::WallBoundarySystem{<:BoundaryModelDummyParticles}, - neighbor::AbstractFluidSystem, u_system, u_neighbor, semi) - # Depending on the density calculator of the boundary model, this NHS is used for - # - kernel summation (`SummationDensity`) - # - continuity equation (`ContinuityDensity`) - # - pressure extrapolation (`AdamiPressureExtrapolation`) - # - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - # The current coordinates of fluids and structured change over time. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], true), - eachindex_y=each_active_particle(neighbor)) -end - -# This function is the same as the one above to avoid ambiguous dispatch when using `Union` -function update_nhs!(neighborhood_search, - system::WallBoundarySystem{<:BoundaryModelDummyParticles}, - neighbor::OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}, - u_system, u_neighbor, semi) - # Depending on the density calculator of the boundary model, this NHS is used for - # - kernel summation (`SummationDensity`) - # - continuity equation (`ContinuityDensity`) - # - pressure extrapolation (`AdamiPressureExtrapolation`) - # - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - # The current coordinates of open boundaries change over time. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], true), - eachindex_y=each_active_particle(neighbor)) -end - -# This function is the same as the one above to avoid ambiguous dispatch when using `Union` -function update_nhs!(neighborhood_search, - system::WallBoundarySystem{<:BoundaryModelDummyParticles}, - neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi) - # Depending on the density calculator of the boundary model, this NHS is used for - # - kernel summation (`SummationDensity`) - # - continuity equation (`ContinuityDensity`) - # - pressure extrapolation (`AdamiPressureExtrapolation`) - # - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - # The current coordinates of fluids and structured change over time. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], true)) -end - -function update_nhs!(neighborhood_search, - system::WallBoundarySystem{<:BoundaryModelDummyParticles}, - neighbor::WallBoundarySystem, - u_system, u_neighbor, semi) - # `system` coordinates only change over time when `system.ismoving[]`. - # `neighbor` coordinates only change over time when `neighbor.ismoving[]`. - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], neighbor.ismoving[])) -end - -function update_nhs!(neighborhood_search, - system::DEMSystem, neighbor::DEMSystem, - u_system, u_neighbor, semi) - # Both coordinates change over time - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true)) -end - -function update_nhs!(neighborhood_search, - system::DEMSystem, neighbor::BoundaryDEMSystem, - u_system, u_neighbor, semi) - # DEM coordinates change over time, the boundary coordinates don't - update!(neighborhood_search, - current_coordinates(u_system, system), - current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, false)) -end - -function update_nhs!(neighborhood_search, - system::WallBoundarySystem, - neighbor::AbstractFluidSystem, - u_system, u_neighbor, semi) - # Don't update. This NHS is never used. - return neighborhood_search -end - -function update_nhs!(neighborhood_search, - system::BoundaryDEMSystem, - neighbor::Union{DEMSystem, BoundaryDEMSystem}, - u_system, u_neighbor, semi) - # Don't update. This NHS is never used. - return neighborhood_search -end - -function update_nhs!(neighborhood_search, - system::Union{WallBoundarySystem, OpenBoundarySystem}, - neighbor::Union{WallBoundarySystem, OpenBoundarySystem}, - u_system, u_neighbor, semi) - # Don't update. This NHS is never used. - return neighborhood_search -end - -# Forward to PointNeighbors.jl -function update!(neighborhood_search, x, y, semi; points_moving=(true, true), - eachindex_y=axes(y, 2)) - PointNeighbors.update!(neighborhood_search, x, y; points_moving, eachindex_y, - parallelization_backend=semi.parallelization_backend) -end - -function check_update_callback(semi) - foreach_system(semi) do system - # This check will be optimized away if the system does not require the callback - if requires_update_callback(system) && !semi.update_callback_used[] - system_name = system |> typeof |> nameof - throw(ArgumentError("`UpdateCallback` is required for `$system_name`")) - end - end -end - -function check_configuration(systems, - nhs::Union{Nothing, PointNeighbors.AbstractNeighborhoodSearch}) - foreach_system(systems) do system - check_configuration(system, systems, nhs) - end - - check_system_color(systems) -end - -check_configuration(system::AbstractSystem, systems, nhs) = nothing - -function check_system_color(systems) - if any(system isa AbstractFluidSystem && !(system isa ParticlePackingSystem) && - !isnothing(system.surface_tension) - for system in systems) - - # System indices of all systems that are either a fluid or a boundary system - system_ids = findall(system isa Union{AbstractFluidSystem, WallBoundarySystem} - for system in systems) - - if length(system_ids) > 1 && sum(i -> systems[i].cache.color, system_ids) == 0 - throw(ArgumentError("If a surface tension model is used the values of at least one system needs to have a color different than 0.")) - end - end -end - -function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs) - if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension) - foreach_system(systems) do neighbor - if neighbor isa AbstractFluidSystem && - isnothing(fluid_system.surface_tension) && - isnothing(fluid_system.surface_normal_method) - throw(ArgumentError("either none or all fluid systems in a simulation need " * - "to use a surface tension model or a surface normal method.")) - end - end - end -end - -function check_configuration(system::WallBoundarySystem, systems, nhs) - (; boundary_model) = system - - foreach_system(systems) do neighbor - if neighbor isa WeaklyCompressibleSPHSystem && - boundary_model isa BoundaryModelDummyParticles && - isnothing(boundary_model.state_equation) - throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without " * - "setting a `state_equation` for all boundary models")) - end - end -end - -function check_configuration(system::TotalLagrangianSPHSystem, systems, nhs) - (; boundary_model) = system - - foreach_system(systems) do neighbor - if neighbor isa AbstractFluidSystem && boundary_model === nothing - throw(ArgumentError("a boundary model for `TotalLagrangianSPHSystem` must be " * - "specified when simulating a fluid-structure interaction.")) - end - end - - if boundary_model isa BoundaryModelDummyParticles && - boundary_model.density_calculator isa ContinuityDensity - throw(ArgumentError("`BoundaryModelDummyParticles` with density calculator " * - "`ContinuityDensity` is not yet supported for a `TotalLagrangianSPHSystem`")) - end -end - -function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, nhs) - foreach_system(systems) do neighbor - if neighbor isa WeaklyCompressibleSPHSystem - throw(ArgumentError("`ImplicitIncompressibleSPHSystem` cannot be used together with - `WeaklyCompressibleSPHSystem`")) - end - end -end - -function check_configuration(system::OpenBoundarySystem, systems, - neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch) - (; boundary_model, boundary_zones) = system - - # Store index of the fluid system. This is necessary for re-linking - # in case we use Adapt.jl to create a new semidiscretization. - fluid_system_index = findfirst(==(system.fluid_system), systems) - system.fluid_system_index[] = fluid_system_index - - if boundary_model isa BoundaryModelCharacteristicsLastiwka && - any(zone -> isnothing(zone.flow_direction), boundary_zones) - throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " * - "Please specify `InFlow()` and `OutFlow()`.")) - end - - if first(PointNeighbors.requires_update(neighborhood_search)) - throw(ArgumentError("`OpenBoundarySystem` requires a neighborhood search " * - "that does not require an update for the first set of coordinates (e.g. `GridNeighborhoodSearch`). " * - "See the PointNeighbors.jl documentation for more details.")) - end -end - -# After `adapt`, the system type information may change. -# This means that systems linking to other systems still point to old systems. -# Therefore, we have to re-link them based on the stored system index. -set_system_links(system, semi) = system - -function set_system_links(system::OpenBoundarySystem, semi) - fluid_system = semi.systems[system.fluid_system_index[]] - - return OpenBoundarySystem(system.boundary_model, - system.initial_condition, - fluid_system, # link to fluid system - system.fluid_system_index, - system.smoothing_kernel, - system.smoothing_length, - system.mass, - system.volume, - system.boundary_candidates, - system.fluid_candidates, - system.boundary_zone_indices, - system.boundary_zones, - system.buffer, - system.pressure_acceleration_formulation, - system.shifting_technique, - system.cache) -end diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl deleted file mode 100644 index 9f1def3b5..000000000 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ /dev/null @@ -1,641 +0,0 @@ -@doc raw""" - BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, - density_calculator, smoothing_kernel, - smoothing_length; viscosity=nothing, - state_equation=nothing, correction=nothing, - reference_particle_spacing=0.0) - -Boundary model for [`WallBoundarySystem`](@ref). - -# Arguments -- `initial_density`: Vector holding the initial density of each boundary particle. -- `hydrodynamic_mass`: Vector holding the "hydrodynamic mass" of each boundary particle. - See description above for more information. -- `density_calculator`: Strategy to compute the hydrodynamic density of the boundary particles. - See description below for more information. -- `smoothing_kernel`: Smoothing kernel should be the same as for the adjacent fluid system. -- `smoothing_length`: Smoothing length should be the same as for the adjacent fluid system. - -# Keywords -- `state_equation`: This should be the same as for the adjacent fluid system - (see e.g. [`StateEquationCole`](@ref)). -- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). -- `viscosity`: Slip (default) or no-slip condition. See description below for further - information. -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -# Examples -```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1) -# Free-slip condition -boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length) - -# No-slip condition -boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length, - viscosity=ViscosityAdami(nu=1e-6)) - -# output -BoundaryModelDummyParticles(AdamiPressureExtrapolation, ViscosityAdami) -``` -""" -struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C} - pressure :: VECTOR # Vector{ELTYPE} - hydrodynamic_mass :: VECTOR # Vector{ELTYPE} - state_equation :: SE - density_calculator :: DC - smoothing_kernel :: K - smoothing_length :: ELTYPE - viscosity :: V - correction :: COR - cache :: C -end - -# The default constructor needs to be accessible for Adapt.jl to work with this struct. -# See the comments in general/gpu.jl for more details. -function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, - density_calculator, smoothing_kernel, - smoothing_length; viscosity=nothing, - state_equation=nothing, correction=nothing, - reference_particle_spacing=0.0) - pressure = initial_boundary_pressure(initial_density, density_calculator, - state_equation) - NDIMS = ndims(smoothing_kernel) - ELTYPE = eltype(smoothing_length) - n_particles = length(initial_density) - - cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., - create_cache_model(initial_density, density_calculator)..., - create_cache_model(correction, initial_density, NDIMS, n_particles)...) - - # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` - if reference_particle_spacing > 0 - # `reference_particle_spacing` has to be set for surface normals to be determined - cache = (; - cache..., # Existing cache fields - reference_particle_spacing=reference_particle_spacing, - initial_colorfield=zeros(ELTYPE, n_particles), - colorfield=zeros(ELTYPE, n_particles), - neighbor_count=zeros(ELTYPE, n_particles)) - end - - return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation, - density_calculator, smoothing_kernel, - smoothing_length, viscosity, correction, cache) -end - -@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel) - -@doc raw""" - AdamiPressureExtrapolation(; pressure_offset=0, allow_loop_flipping=true) - -`density_calculator` for `BoundaryModelDummyParticles`. - -# Keywords -- `pressure_offset=0`: Sometimes it is necessary to artificially increase the boundary pressure - to prevent penetration, which is possible by increasing this value. -- `allow_loop_flipping=true`: Allow to flip the loop order for the pressure extrapolation. - Disable to prevent error variations between simulations with - different numbers of threads. - Usually, the first (multithreaded) loop is over the boundary - particles and the second loop over the fluid neighbors. - When the number of boundary particles is larger than - `ceil(0.5 * nthreads())` times the number of fluid particles, - it is usually more efficient to flip the loop order and loop - over the fluid particles first. - The factor depends on the number of threads, as the flipped - loop is not thread parallelizable. - This can cause error variations between simulations with - different numbers of threads. - -""" -struct AdamiPressureExtrapolation{ELTYPE} - pressure_offset :: ELTYPE - allow_loop_flipping :: Bool - - function AdamiPressureExtrapolation(; pressure_offset=0, allow_loop_flipping=true) - return new{eltype(pressure_offset)}(pressure_offset, allow_loop_flipping) - end -end - -@doc raw""" - BernoulliPressureExtrapolation(; pressure_offset=0, factor=1) - -`density_calculator` for `BoundaryModelDummyParticles`. - -# Keywords -- `pressure_offset=0`: Sometimes it is necessary to artificially increase the boundary pressure - to prevent penetration, which is possible by increasing this value. -- `factor=1` : Setting `factor` allows to just increase the strength of the dynamic - pressure part. -- `allow_loop_flipping=true`: Allow to flip the loop order for the pressure extrapolation. - Disable to prevent error variations between simulations with - different numbers of threads. - Usually, the first (multithreaded) loop is over the boundary - particles and the second loop over the fluid neighbors. - When the number of boundary particles is larger than - `ceil(0.5 * nthreads())` times the number of fluid particles, - it is usually more efficient to flip the loop order and loop - over the fluid particles first. - The factor depends on the number of threads, as the flipped - loop is not thread parallelizable. - This can cause error variations between simulations with - different numbers of threads. - -""" -struct BernoulliPressureExtrapolation{ELTYPE} - pressure_offset :: ELTYPE - factor :: ELTYPE - allow_loop_flipping :: Bool - - function BernoulliPressureExtrapolation(; pressure_offset=0, factor=1, - allow_loop_flipping=true) - return new{eltype(pressure_offset)}(pressure_offset, factor, allow_loop_flipping) - end -end - -@doc raw""" - PressureMirroring() - -`density_calculator` for `BoundaryModelDummyParticles`. - -!!! note - This boundary model requires high viscosity for stability with WCSPH. - It also produces significantly worse results than [`AdamiPressureExtrapolation`](@ref) - and is not more efficient because smaller time steps are required due to more noise - in the pressure. - We added this model only for research purposes and for comparison with - [SPlisHSPlasH](https://github.com/InteractiveComputerGraphics/SPlisHSPlasH). -""" -struct PressureMirroring end - -@doc raw""" - PressureZeroing() - -`density_calculator` for `BoundaryModelDummyParticles`. - -!!! note - This boundary model produces significantly worse results than all other models and - is only included for research purposes. -""" -struct PressureZeroing end - -@inline create_cache_model(correction, density, NDIMS, nparticles) = (;) - -function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particles) - return (; kernel_correction_coefficient=similar(density)) -end - -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 - -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 - -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, - ::Union{SummationDensity, PressureMirroring, PressureZeroing}) - density = copy(initial_density) - - return (; density) -end - -@inline create_cache_model(initial_density, ::ContinuityDensity) = (; initial_density) - -function create_cache_model(initial_density, - ::Union{AdamiPressureExtrapolation, - BernoulliPressureExtrapolation}) - density = copy(initial_density) - volume = similar(initial_density) - - return (; density, volume) -end - -@inline create_cache_model(viscosity::Nothing, n_particles, n_dims) = (;) - -function create_cache_model(viscosity, n_particles, n_dims) - ELTYPE = eltype(viscosity.epsilon) - - wall_velocity = zeros(ELTYPE, n_dims, n_particles) - - return (; wall_velocity) -end - -@inline reset_cache!(cache, viscosity) = set_zero!(cache.volume) - -function reset_cache!(cache, viscosity::ViscosityAdami) - (; volume, wall_velocity) = cache - - set_zero!(volume) - set_zero!(wall_velocity) - - 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 current_density(v, model::BoundaryModelDummyParticles, system) - return current_density(v, model.density_calculator, model) -end - -@inline function current_density(v, - ::Union{SummationDensity, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, - BernoulliPressureExtrapolation}, - model::BoundaryModelDummyParticles) - # When using `SummationDensity`, the density is stored in the cache - return model.cache.density -end - -@inline function current_density(v, ::ContinuityDensity, - model::BoundaryModelDummyParticles) - # When using `ContinuityDensity`, the density is stored in the last row of `v` - return view(v, size(v, 1), :) -end - -@inline function current_pressure(v, model::BoundaryModelDummyParticles, system) - return model.pressure -end - -@inline function update_density!(boundary_model::BoundaryModelDummyParticles, - system, v, u, v_ode, u_ode, semi) - (; density_calculator) = boundary_model - - compute_density!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) - - return boundary_model -end - -function compute_density!(boundary_model, - ::Union{ContinuityDensity, AdamiPressureExtrapolation, - BernoulliPressureExtrapolation, - PressureMirroring, PressureZeroing}, - system, v, u, v_ode, u_ode, semi) - # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`. - # For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. - return boundary_model -end - -@inline function update_pressure!(boundary_model::BoundaryModelDummyParticles, - system, v, u, v_ode, u_ode, semi) - (; correction, density_calculator) = boundary_model - - compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) - - # These are only computed when using corrections - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) - compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode, - semi) - # `kernel_correct_density!` only performed for `SummationDensity` - kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction, - density_calculator) - - return boundary_model -end - -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) = boundary_model - (; correction_matrix) = cache - - system_coords = current_coordinates(u, system) - - compute_gradient_correction_matrix!(correction_matrix, system, system_coords, - v_ode, u_ode, semi, correction, smoothing_kernel) -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 - - summation_density!(system, semi, u, u_ode, density, particles=eachparticle(system)) -end - -function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityDensity}, - system, v, u, v_ode, u_ode, semi) - - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). - @threaded semi for particle in eachparticle(system) - apply_state_equation!(boundary_model, current_density(v, system, 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) -end - -function compute_pressure!(boundary_model, - ::Union{AdamiPressureExtrapolation, - BernoulliPressureExtrapolation}, - system, v, u, v_ode, u_ode, semi) - (; pressure, cache, viscosity) = boundary_model - (; allow_loop_flipping) = boundary_model.density_calculator - - set_zero!(pressure) - - # Set `volume` to zero. For `ViscosityAdami` the `wall_velocity` is also set to zero. - reset_cache!(cache, viscosity) - - system_coords = current_coordinates(u, system) - - # Use all other systems for the pressure extrapolation - @trixi_timeit timer() "compute boundary pressure" foreach_system(semi) do neighbor_system - v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) - u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - # This is an optimization for simulations with large and complex boundaries. - # Especially in 3D simulations with large and/or complex structures outside - # of areas with permanent flow. - # Note: The version iterating neighbors first is not thread-parallelizable - # and thus not GPU-compatible. - # The factor is based on the achievable speed-up of the thread parallelizable version. - # Use the parallel version if the number of boundary particles is not much larger - # than the number of fluid particles. - n_boundary_particles = nparticles(system) - n_fluid_particles = nparticles(neighbor_system) - speedup = ceil(Int, Threads.nthreads() / 2) - is_gpu = system_coords isa AbstractGPUArray - condition_boundary = n_boundary_particles < speedup * n_fluid_particles - parallelize = is_gpu || condition_boundary || !allow_loop_flipping - - # Loop over boundary particles and then the neighboring fluid particles - # to extrapolate fluid pressure to the boundaries. - boundary_pressure_extrapolation!(Val(parallelize), boundary_model, system, - neighbor_system, system_coords, neighbor_coords, v, - v_neighbor_system, semi) - end - - @trixi_timeit timer() "inverse state equation" @threaded semi for particle in - eachparticle(system) - compute_adami_density!(boundary_model, system, v, 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, v, 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 @inbounds volume[particle] > eps() - @inbounds pressure[particle] /= volume[particle] - - # To impose no-slip condition - compute_wall_velocity!(viscosity, system, v, particle) - end - - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). - @inbounds pressure[particle] = max(pressure[particle], 0) - - # 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}, - system, v, u, v_ode, u_ode, semi) - # No pressure update needed with `PressureMirroring` and `PressureZeroing`. - return boundary_model -end - -@inline function boundary_pressure_extrapolation!(parallel, boundary_model, system, - neighbor_system, system_coords, - neighbor_coords, v, v_neighbor_system, - semi) - return boundary_model -end - -@inline function boundary_pressure_extrapolation!(parallel::Val{true}, boundary_model, - system, - neighbor_system::Union{AbstractFluidSystem, - OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - system_coords, neighbor_coords, v, - v_neighbor_system, semi) - (; pressure, cache, viscosity, density_calculator) = boundary_model - (; pressure_offset) = density_calculator - - # Loop over all pairs of particles and neighbors within the kernel cutoff - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; - points=eachparticle(system)) do particle, neighbor, - pos_diff, distance - boundary_pressure_inner!(boundary_model, density_calculator, system, - neighbor_system, v, v_neighbor_system, particle, neighbor, - pos_diff, distance, viscosity, cache, pressure, - pressure_offset) - end -end - -# Loop over fluid particles and then the neighboring boundary particles -# to extrapolate fluid pressure to the boundaries. -# Note that this needs to be serial, as we are writing into the same -# pressure entry from different loop iterations. -@inline function boundary_pressure_extrapolation!(parallel::Val{false}, boundary_model, - system, - neighbor_system::Union{AbstractFluidSystem, - OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - system_coords, neighbor_coords, - v, v_neighbor_system, semi) - (; pressure, cache, viscosity, density_calculator) = boundary_model - (; pressure_offset) = density_calculator - - # This needs to be serial to avoid race conditions when writing into `system` - foreach_point_neighbor(neighbor_system, system, neighbor_coords, system_coords, semi; - points=each_integrated_particle(neighbor_system), - parallelization_backend=SerialBackend()) do neighbor, particle, - pos_diff, distance - # Since neighbor and particle are switched - pos_diff = -pos_diff - boundary_pressure_inner!(boundary_model, density_calculator, system, - neighbor_system, v, v_neighbor_system, particle, neighbor, - pos_diff, distance, viscosity, cache, pressure, - pressure_offset) - end -end - -@inline function boundary_pressure_inner!(boundary_model, boundary_density_calculator, - system, - neighbor_system::Union{AbstractFluidSystem, - OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - v, v_neighbor_system, particle, neighbor, - pos_diff, distance, viscosity, cache, pressure, - pressure_offset) - density_neighbor = @inbounds current_density(v_neighbor_system, neighbor_system, - neighbor) - - # Fluid pressure term - fluid_pressure = @inbounds current_pressure(v_neighbor_system, neighbor_system, - neighbor) - - # Hydrostatic pressure term from fluid and boundary acceleration - # TODO: rename acceleration and add a function `system_external_acceleration` - resulting_acceleration = @inbounds acceleration_source(neighbor_system) - - @inbounds current_acceleration(system, particle) - hydrostatic_pressure = dot(resulting_acceleration, density_neighbor * pos_diff) - - # Additional dynamic pressure term (only with `BernoulliPressureExtrapolation`) - dynamic_pressure_ = dynamic_pressure(boundary_density_calculator, density_neighbor, - v, v_neighbor_system, pos_diff, distance, - particle, neighbor, system, neighbor_system) - - sum_pressures = pressure_offset + fluid_pressure + dynamic_pressure_ + - hydrostatic_pressure - - kernel_weight = smoothing_kernel(boundary_model, distance, particle) - - @inbounds pressure[particle] += sum_pressures * kernel_weight - @inbounds cache.volume[particle] += kernel_weight - - compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, - kernel_weight, particle, neighbor) -end - -@inline function dynamic_pressure(boundary_density_calculator, density_neighbor, v, - v_neighbor_system, pos_diff, distance, particle, neighbor, - system, neighbor_system) - return zero(density_neighbor) -end - -@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, - density_neighbor, v, v_neighbor_system, pos_diff, - distance, particle, neighbor, - system::AbstractBoundarySystem, neighbor_system) - if system.ismoving[] - relative_velocity = current_velocity(v, system, particle) .- - current_velocity(v_neighbor_system, neighbor_system, neighbor) - normal_velocity = dot(relative_velocity, pos_diff) - - return boundary_density_calculator.factor * density_neighbor * - normal_velocity^2 / distance / 2 - end - return zero(density_neighbor) -end - -@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, - density_neighbor, v, v_neighbor_system, pos_diff, - distance, particle, neighbor, - system::AbstractStructureSystem, neighbor_system) - relative_velocity = current_velocity(v, system, particle) .- - current_velocity(v_neighbor_system, neighbor_system, neighbor) - normal_velocity = dot(relative_velocity, pos_diff) / distance - - return boundary_density_calculator.factor * density_neighbor * - dot(normal_velocity, normal_velocity) / 2 -end - -function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, - kernel_weight, particle, neighbor) - return cache -end - -function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, neighbor_system, - v_neighbor_system, kernel_weight, particle, neighbor) - v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - - for dim in eachindex(v_b) - @inbounds cache.wall_velocity[dim, particle] += kernel_weight * v_b[dim] - end - - return cache -end - -@inline function compute_wall_velocity!(viscosity::Nothing, system, system_coords, particle) - return viscosity -end - -@inline function compute_wall_velocity!(viscosity, system, v, particle) - (; boundary_model) = system - (; cache) = boundary_model - (; volume, wall_velocity) = cache - - # Prescribed velocity of the boundary particle. - # This velocity is zero when not using moving boundaries. - v_boundary = current_velocity(v, system, particle) - - for dim in eachindex(v_boundary) - # The second term is the precalculated smoothed velocity field of the fluid - new_velocity = @inbounds 2 * v_boundary[dim] - - wall_velocity[dim, particle] / volume[particle] - @inbounds wall_velocity[dim, particle] = new_velocity - end - return viscosity -end - -@inline function inverse_state_equation!(density, state_equation, pressure, particle) - @inbounds density[particle] = inverse_state_equation(state_equation, pressure[particle]) - return density -end - -@inline function inverse_state_equation!(density, state_equation::Nothing, pressure, - particle) - # The density is constant when using EDAC - return density -end - -@inline function correction_matrix(system::AbstractBoundarySystem, particle) - extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) -end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl deleted file mode 100644 index cea18484d..000000000 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ /dev/null @@ -1,134 +0,0 @@ -# Fluid-fluid and fluid-boundary interaction -function interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - particle_system::EntropicallyDampedSPHSystem, - neighbor_system, semi) - (; sound_speed, density_calculator, correction, nu_edac) = particle_system - - system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - surface_tension_a = surface_tension_model(particle_system) - surface_tension_b = surface_tension_model(neighbor_system) - - # Loop over all pairs of particles and neighbors within the kernel cutoff - foreach_point_neighbor(particle_system, neighbor_system, - system_coords, neighbor_coords, semi; - points=each_integrated_particle(particle_system)) do particle, - neighbor, - pos_diff, - distance - # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are - # in bounds of the respective system. For performance reasons, we use `@inbounds` - # in this hot loop to avoid bounds checking when extracting particle quantities. - rho_a = @inbounds current_density(v_particle_system, particle_system, particle) - rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor) - - p_a = @inbounds current_pressure(v_particle_system, particle_system, particle) - p_b = @inbounds current_pressure(v_neighbor_system, neighbor_system, neighbor) - - # This technique by Basa et al. 2017 (10.1002/fld.1927) aims to reduce numerical - # errors due to large pressures by subtracting the average pressure of neighboring - # particles. - # It results in significant improvement for EDAC, especially with TVF, - # but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2. - # Note that the return value is zero when not using average pressure reduction. - p_avg = @inbounds average_pressure(particle_system, particle) - - m_a = @inbounds hydrodynamic_mass(particle_system, particle) - m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) - - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - - dv_pressure = pressure_acceleration(particle_system, neighbor_system, - particle, neighbor, - m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, - rho_b, pos_diff, distance, grad_kernel, - correction) - - dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) - - # Extra terms in the momentum equation when using a shifting technique - dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), - particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, m_a, m_b, rho_a, rho_b, - pos_diff, distance, grad_kernel, correction) - - dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) - - dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, - particle, neighbor, pos_diff, distance) - - dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension + - dv_adhesion - - for i in 1:ndims(particle_system) - @inbounds dv[i, particle] += dv_particle[i] - end - - v_diff = current_velocity(v_particle_system, particle_system, particle) - - current_velocity(v_neighbor_system, neighbor_system, neighbor) - - pressure_evolution!(dv, particle_system, neighbor_system, v_diff, grad_kernel, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, p_a, p_b, rho_a, rho_b, nu_edac) - - # TODO If variable smoothing_length is used, this should use the neighbor smoothing length - # Propagate `@inbounds` to the continuity equation, which accesses particle data - @inbounds continuity_equation!(dv, density_calculator, particle_system, - neighbor_system, v_particle_system, - v_neighbor_system, particle, neighbor, - pos_diff, distance, m_b, rho_a, rho_b, grad_kernel) - end - - return dv -end - -@inline function pressure_evolution!(dv, particle_system, neighbor_system, v_diff, - grad_kernel, particle, neighbor, - pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b, nu_edac) - volume_a = m_a / rho_a - volume_b = m_b / rho_b - volume_term = (volume_a^2 + volume_b^2) / m_a - - # EDAC pressure evolution - pressure_diff = p_a - p_b - - # This is basically the continuity equation times `sound_speed^2` - artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) - - eta_a = rho_a * nu_edac - eta_b = rho_b * nu_edac - eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) - - smoothing_length_average = (smoothing_length(particle_system, particle) + - smoothing_length(neighbor_system, neighbor)) / 2 - tmp = eta_tilde / (distance^2 + smoothing_length_average^2 / 100) - - # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 - # They argued that the formulation is more flexible because of the possibility to formulate - # different inter-particle averages or to assume different inter-particle distributions. - # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. - # - # TODO: Is there a better formulation to discretize the Laplace operator? - # Because when using this formulation for the pressure acceleration, it is not - # energy conserving. - # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 - # - # This is similar to density diffusion in WCSPH - damping_term = volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) - - # Pressure is stored in `v` right after the velocity - dv[ndims(particle_system) + 1, particle] += artificial_eos + damping_term - - return dv -end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl deleted file mode 100644 index 9368e7f71..000000000 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ /dev/null @@ -1,382 +0,0 @@ -@doc raw""" - EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed; - pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), - shifting_technique=nothing, - alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, - reference_particle_spacing=0.0, source_terms=nothing) - -System for particles of a fluid. -As opposed to the [weakly compressible SPH scheme](@ref wcsph), which uses an equation of state, -this scheme uses a pressure evolution equation to calculate the pressure. -See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more details on the method. - -# Arguments -- `initial_condition`: Initial condition representing the system's particles. -- `sound_speed`: Speed of sound. -- `smoothing_kernel`: Smoothing kernel to be used for this system. - See [Smoothing Kernels](@ref smoothing_kernel). -- `smoothing_length`: Smoothing length to be used for this system. - See [Smoothing Kernels](@ref smoothing_kernel). - -# Keywords -- `viscosity`: Viscosity model for this system (default: no viscosity). - Recommended: [`ViscosityAdami`](@ref). -- `acceleration`: Acceleration vector for the system. (default: zero vector) -- `pressure_acceleration`: Pressure acceleration formulation (default: inter-particle averaged pressure). - When set to `nothing`, the pressure acceleration formulation for the - corresponding [density calculator](@ref density_calculator) is chosen. -- `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) -- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity - formulation](@ref transport_velocity_formulation) to use - with this system. Default is no shifting. -- `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles - from the local pressure (default: `true` when using shifting, `false` otherwise). -- `buffer_size`: Number of buffer particles. - This is needed when simulating with [`OpenBoundarySystem`](@ref). -- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) -- `source_terms`: Additional source terms for this system. Has to be either `nothing` - (by default), or a function of `(coords, velocity, density, pressure, t)` - (which are the quantities of a single particle), returning a `Tuple` - or `SVector` that is to be added to the acceleration of that particle. - See, for example, [`SourceTermDamping`](@ref). - Note that these source terms will not be used in the calculation of the - boundary pressure when using a boundary with - [`BoundaryModelDummyParticles`](@ref) and [`AdamiPressureExtrapolation`](@ref). - The keyword argument `acceleration` should be used instead for - gravity-like source terms. -- `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) -- `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: The value used to initialize the color of particles in the system. - -""" -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV, - AVGP, ST, SRFT, SRFN, B, PR, - C} <: AbstractFluidSystem{NDIMS} - initial_condition :: IC - mass :: M # Vector{ELTYPE}: [particle] - density_calculator :: DC - smoothing_kernel :: K - sound_speed :: ELTYPE - viscosity :: V - nu_edac :: ELTYPE - acceleration :: SVector{NDIMS, ELTYPE} - correction :: COR - pressure_acceleration_formulation :: PF - shifting_technique :: TV - average_pressure_reduction :: AVGP - source_terms :: ST - surface_tension :: SRFT - surface_normal_method :: SRFN - buffer :: B - particle_refinement :: PR - cache :: C -end - -# The default constructor needs to be accessible for Adapt.jl to work with this struct. -# See the comments in general/gpu.jl for more details. -function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed; - pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), - shifting_technique=nothing, - average_pressure_reduction=(!isnothing(shifting_technique)), - alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel)), - correction=nothing, - source_terms=nothing, surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, - reference_particle_spacing=0.0, color_value=1) - buffer = isnothing(buffer_size) ? nothing : - SystemBuffer(nparticles(initial_condition), buffer_size) - - particle_refinement = nothing # TODO - - initial_condition = allocate_buffer(initial_condition, buffer) - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - mass = copy(initial_condition.mass) - n_particles = length(initial_condition.mass) - - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end - - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end - - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end - - if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) - end - - if correction isa ShepardKernelCorrection && - density_calculator isa ContinuityDensity - throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) - end - - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - correction) - - avg_pressure_reduction = Val(average_pressure_reduction) - - nu_edac = (alpha * smoothing_length * sound_speed) / 8 - - cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_shifting(initial_condition, shifting_technique)..., - create_cache_avg_pressure_reduction(initial_condition, - avg_pressure_reduction)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., - create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, - n_particles)..., - create_cache_refinement(initial_condition, particle_refinement, - smoothing_length)..., - create_cache_correction(correction, initial_condition.density, NDIMS, - n_particles)..., - color=Int(color_value)) - - # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` - if reference_particle_spacing > 0 - # `reference_particle_spacing` has to be set for surface normals to be determined - cache = (; - cache..., # Existing cache fields - reference_particle_spacing=reference_particle_spacing) - end - - EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), - typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), typeof(shifting_technique), - typeof(avg_pressure_reduction), typeof(source_terms), - typeof(surface_tension), typeof(surface_normal_method), - typeof(buffer), Nothing, - typeof(cache)}(initial_condition, mass, density_calculator, - smoothing_kernel, sound_speed, viscosity, - nu_edac, acceleration_, correction, - pressure_acceleration, shifting_technique, - avg_pressure_reduction, - source_terms, surface_tension, - surface_normal_method, buffer, - particle_refinement, cache) -end - -create_cache_avg_pressure_reduction(initial_condition, ::Val{false}) = (;) - -function create_cache_avg_pressure_reduction(initial_condition, ::Val{true}) - pressure_average = copy(initial_condition.pressure) - neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) - - return (; pressure_average, neighbor_counter) -end - -function Base.show(io::IO, system::EntropicallyDampedSPHSystem) - @nospecialize system # reduce precompilation time - - print(io, "EntropicallyDampedSPHSystem{", ndims(system), "}(") - print(io, system.density_calculator) - print(io, ", ", system.correction) - print(io, ", ", system.viscosity) - print(io, ", ", system.smoothing_kernel) - print(io, ", ", system.acceleration) - print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) - print(io, ") with ", nparticles(system), " particles") -end - -function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSystem) - @nospecialize system # reduce precompilation time - - if get(io, :compact, false) - show(io, system) - else - summary_header(io, "EntropicallyDampedSPHSystem{$(ndims(system))}") - if system.buffer isa SystemBuffer - summary_line(io, "#particles", nparticles(system)) - summary_line(io, "#buffer_particles", system.buffer.buffer_size) - else - summary_line(io, "#particles", nparticles(system)) - end - summary_line(io, "density calculator", - system.density_calculator |> typeof |> nameof) - summary_line(io, "correction method", - system.correction |> typeof |> nameof) - summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) - summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") - summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "shifting technique", system.shifting_technique) - summary_line(io, "average pressure reduction", - typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") - summary_line(io, "acceleration", system.acceleration) - summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) - summary_footer(io) - end -end - -@inline function Base.eltype(::EntropicallyDampedSPHSystem{<:Any, ELTYPE}) where {ELTYPE} - return ELTYPE -end - -@inline function v_nvariables(system::EntropicallyDampedSPHSystem) - return v_nvariables(system, system.density_calculator) -end - -@inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::SummationDensity) - return ndims(system) + 1 -end - -@inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::ContinuityDensity) - return ndims(system) + 2 -end - -@inline buffer(system::EntropicallyDampedSPHSystem) = system.buffer - -system_correction(system::EntropicallyDampedSPHSystem) = system.correction - -@inline function current_velocity(v, system::EntropicallyDampedSPHSystem) - return view(v, 1:ndims(system), :) -end - -@inline system_state_equation(system::EntropicallyDampedSPHSystem) = nothing - -@inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed - -@inline shifting_technique(system::EntropicallyDampedSPHSystem) = system.shifting_technique - -@propagate_inbounds function average_pressure(system::EntropicallyDampedSPHSystem, particle) - average_pressure(system, system.average_pressure_reduction, particle) -end - -@propagate_inbounds function average_pressure(system, ::Val{true}, particle) - return system.cache.pressure_average[particle] -end - -@inline average_pressure(system, ::Val{false}, particle) = zero(eltype(system)) - -@inline function current_density(v, system::EntropicallyDampedSPHSystem) - return current_density(v, system.density_calculator, system) -end - -@inline function current_density(v, ::SummationDensity, - system::EntropicallyDampedSPHSystem) - # When using `SummationDensity`, the density is stored in the cache - return system.cache.density -end - -@inline function current_density(v, ::ContinuityDensity, - system::EntropicallyDampedSPHSystem) - # When using `ContinuityDensity`, the density is stored in the last row of `v` - return view(v, size(v, 1), :) -end - -@inline function current_pressure(v, system::EntropicallyDampedSPHSystem) - return view(v, ndims(system) + 1, :) -end - -function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, - v_ode, u_ode, semi, t) - compute_density!(system, u, u_ode, semi, system.density_calculator) -end - -function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) - compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, - t) - compute_surface_delta_function!(system, system.surface_tension, semi) -end - -function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) - (; surface_tension) = system - - # Surface normal of neighbor and boundary needs to have been calculated already - compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi) - update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) -end - -# No average pressure reduction is used -function update_average_pressure!(system, ::Val{false}, v_ode, u_ode, semi) - return system -end - -# This technique by Basa et al. 2017 (10.1002/fld.1927) aims to reduce numerical errors -# due to large pressures by subtracting the average pressure of neighboring particles. -# It results in significant improvement for EDAC, especially with TVF, -# but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2. -# See eq. 16 and 17 in Ramachandran & Puri (2019) for an explanation of the technique. -function update_average_pressure!(system, ::Val{true}, v_ode, u_ode, semi) - (; cache) = system - (; pressure_average, neighbor_counter) = cache - - set_zero!(pressure_average) - set_zero!(neighbor_counter) - - u = wrap_u(u_ode, system, semi) - - # Use all other systems for the average pressure - @trixi_timeit timer() "compute average pressure" 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) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - # Loop over all pairs of particles and neighbors within the kernel cutoff. - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_integrated_particle(system)) do particle, - neighbor, - pos_diff, - distance - pressure_average[particle] += current_pressure(v_neighbor_system, - neighbor_system, neighbor) - neighbor_counter[particle] += 1 - end - end - - # We do not need to check for zero division here, as `neighbor_counter = 1` - # for zero neighbors. That is, the `particle` itself is also taken into account. - pressure_average ./= neighbor_counter - - return system -end - -function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::SummationDensity) - # Note that `.=` is very slightly faster, but not GPU-compatible - v0[ndims(system) + 1, :] = system.initial_condition.pressure - - return v0 -end - -function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::ContinuityDensity) - # Note that `.=` is very slightly faster, but not GPU-compatible - v0[end, :] = system.initial_condition.density - v0[ndims(system) + 1, :] = system.initial_condition.pressure - - return v0 -end - -function restart_with!(system::EntropicallyDampedSPHSystem, v, u) - for particle in each_integrated_particle(system) - system.initial_condition.coordinates[:, particle] .= u[:, particle] - system.initial_condition.velocity[:, particle] .= v[1:ndims(system), particle] - system.initial_condition.pressure[particle] = v[end, particle] - end -end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl deleted file mode 100644 index 5ca225970..000000000 --- a/src/schemes/fluid/fluid.jl +++ /dev/null @@ -1,249 +0,0 @@ -# WARNING! -# These functions are intended to be used internally to set the density -# of newly activated particles in a callback. -# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` -# outside of callbacks. -@inline function set_particle_density!(v, system::AbstractFluidSystem, particle, density) - current_density(v, system)[particle] = density - - return v -end - -# WARNING! -# These functions are intended to be used internally to set the pressure -# of newly activated particles in a callback. -# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` -# outside of callbacks. -@inline function set_particle_pressure!(v, system::AbstractFluidSystem, particle, pressure) - current_pressure(v, system)[particle] = pressure - - return v -end - -function create_cache_density(initial_condition, ::SummationDensity) - density = similar(initial_condition.density) - - return (; density) -end - -function create_cache_density(ic, ::ContinuityDensity) - # Density in this case is added to the end of `v` and allocated by modifying `v_nvariables`. - return (;) -end - -function create_cache_refinement(initial_condition, ::Nothing, smoothing_length) - smoothing_length_factor = smoothing_length / initial_condition.particle_spacing - return (; smoothing_length, smoothing_length_factor) -end - -# TODO -function create_cache_refinement(initial_condition, refinement, smoothing_length) - # TODO: If refinement is not `Nothing` and `correction` is not `Nothing`, then throw an error -end - -@propagate_inbounds function hydrodynamic_mass(system::AbstractFluidSystem, particle) - return system.mass[particle] -end - -function smoothing_length(system::AbstractFluidSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::AbstractFluidSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function initial_smoothing_length(system::AbstractFluidSystem) - return initial_smoothing_length(system, system.particle_refinement) -end - -initial_smoothing_length(system, ::Nothing) = system.cache.smoothing_length - -function initial_smoothing_length(system, refinement) - # TODO - return system.cache.initial_smoothing_length_factor * - system.initial_condition.particle_spacing -end - -@inline function particle_spacing(system::AbstractFluidSystem, particle) - return particle_spacing(system, system.particle_refinement, particle) -end - -@inline particle_spacing(system, ::Nothing, _) = system.initial_condition.particle_spacing - -@inline function particle_spacing(system, refinement, particle) - (; smoothing_length_factor) = system.cache - return smoothing_length(system, particle) / smoothing_length_factor -end - -function write_u0!(u0, system::AbstractFluidSystem) - (; initial_condition) = system - - # This is as fast as a loop with `@inbounds`, but it's GPU-compatible - indices = CartesianIndices(initial_condition.coordinates) - copyto!(u0, indices, initial_condition.coordinates, indices) - - return u0 -end - -function write_v0!(v0, system::AbstractFluidSystem) - # This is as fast as a loop with `@inbounds`, but it's GPU-compatible - indices = CartesianIndices(system.initial_condition.velocity) - copyto!(v0, indices, system.initial_condition.velocity, indices) - - write_v0!(v0, system, system.density_calculator) - - return v0 -end - -write_v0!(v0, system::AbstractFluidSystem, _) = v0 - -# To account for boundary effects in the viscosity term of the RHS, use the viscosity model -# of the neighboring particle systems. - -@inline function viscosity_model(system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem) - return neighbor_system.viscosity -end - -@inline function viscosity_model(system::AbstractFluidSystem, - neighbor_system::AbstractBoundarySystem) - return neighbor_system.boundary_model.viscosity -end - -@inline system_state_equation(system::AbstractFluidSystem) = system.state_equation - -@inline acceleration_source(system::AbstractFluidSystem) = system.acceleration - -function compute_density!(system, u, u_ode, semi, ::ContinuityDensity) - # No density update with `ContinuityDensity` - return system -end - -function compute_density!(system, u, u_ode, semi, ::SummationDensity) - (; cache) = system - (; density) = cache # Density is in the cache for SummationDensity - - summation_density!(system, semi, u, u_ode, density) -end - -# With 'SummationDensity', density is calculated in wcsph/system.jl:compute_density! -@inline function continuity_equation!(dv, density_calculator::SummationDensity, - particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, grad_kernel) - return dv -end - -# This formulation was chosen to be consistent with the used pressure_acceleration formulations -@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity, - particle_system::AbstractFluidSystem, - neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, grad_kernel) - vdiff = current_velocity(v_particle_system, particle_system, particle) - - current_velocity(v_neighbor_system, neighbor_system, neighbor) - - dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel) - - # Artificial density diffusion should only be applied to systems representing a fluid - # with the same physical properties i.e. density and viscosity. - # TODO: shouldn't be applied to particles on the interface (depends on PR #539) - if particle_system === neighbor_system - density_diffusion!(dv, density_diffusion(particle_system), - v_particle_system, particle, neighbor, - pos_diff, distance, m_b, rho_a, rho_b, particle_system, - grad_kernel) - end - - continuity_equation_shifting!(dv, shifting_technique(particle_system), - particle_system, neighbor_system, - particle, neighbor, grad_kernel, rho_a, rho_b, m_b) -end - -function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, semi) - (; viscosity, acceleration, surface_tension) = system - - # TODO - smoothing_length_ = initial_smoothing_length(system) - - dt_viscosity = 0.125 * smoothing_length_^2 - if !isnothing(system.viscosity) - dt_viscosity = dt_viscosity / - kinematic_viscosity(system, viscosity, smoothing_length_, - system_sound_speed(system)) - end - - # TODO Adami et al. (2012) just use the gravity here, but Antuono et al. (2012) - # are using a per-particle acceleration. Is that supposed to be the previous RHS? - # Morris (2000) also uses the acceleration and cites Monaghan (1992) - dt_acceleration = 0.25 * sqrt(smoothing_length_ / norm(acceleration)) - - # TODO Everyone seems to be doing this differently. - # Morris (2000) uses the additional condition CFL < 0.25. - # Sun et al. (2017) only use h / c (because c depends on v_max as c >= 10 v_max). - # Adami et al. (2012) use h / (c + v_max) with a fixed CFL of 0.25. - # Antuono et al. (2012) use h / (c + v_max + h * pi_max), where pi is the viscosity coefficient. - # Antuono et al. (2015) use h / (c + h * pi_max). - # - # See docstring of the callback for the references. - dt_sound_speed = cfl_number * smoothing_length_ / system_sound_speed(system) - - # Eq. 28 in Morris (2000) - dt = min(dt_viscosity, dt_acceleration, dt_sound_speed) - if surface_tension isa SurfaceTensionMorris || - surface_tension isa SurfaceTensionMomentumMorris - v = wrap_v(v_ode, system, semi) - dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 / - (2 * pi * surface_tension.surface_tension_coefficient)) - dt = min(dt, dt_surface_tension) - end - - return dt -end - -@inline function surface_tension_model(system::AbstractFluidSystem) - return system.surface_tension -end - -@inline function surface_tension_model(system) - return nothing -end - -@inline function surface_normal_method(system::AbstractFluidSystem) - return system.surface_normal_method -end - -@inline function surface_normal_method(system) - return nothing -end - -function system_data(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi) - (; mass) = system - - dv = wrap_v(dv_ode, system, semi) - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - coordinates = current_coordinates(u, system) - velocity = current_velocity(v, system) - acceleration = current_velocity(dv, system) - density = current_density(v, system) - pressure = current_pressure(v, system) - - return (; coordinates, velocity, mass, density, pressure, acceleration) -end - -function available_data(::AbstractFluidSystem) - return (:coordinates, :velocity, :mass, :density, :pressure, :acceleration) -end - -include("pressure_acceleration.jl") -include("viscosity.jl") -include("shifting_techniques.jl") -include("surface_tension.jl") -include("surface_normal_sph.jl") -include("weakly_compressible_sph/weakly_compressible_sph.jl") -include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl deleted file mode 100644 index 7a0232c51..000000000 --- a/src/schemes/fluid/surface_tension.jl +++ /dev/null @@ -1,322 +0,0 @@ -abstract type AbstractSurfaceTension end -abstract type AkinciTypeSurfaceTension <: AbstractSurfaceTension end - -@doc raw""" - CohesionForceAkinci(surface_tension_coefficient=1.0) - -This model only implements the cohesion force of the Akinci [Akinci2013](@cite) surface tension model. - -See [`surface_tension`](@ref) for more details. - -# Keywords -- `surface_tension_coefficient=1.0`: Modifies the intensity of the surface tension-induced force, - enabling the tuning of the fluid's surface tension properties within the simulation. -""" -struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension - surface_tension_coefficient::ELTYPE - - function CohesionForceAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) - end -end - -@doc raw""" - SurfaceTensionAkinci(surface_tension_coefficient=1.0) - -Implements a model for surface tension and adhesion effects drawing upon the -principles outlined by Akinci [Akinci2013](@cite). This model is instrumental in capturing the nuanced -behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or -separation, by utilizing intra-particle forces. - -See [`surface_tension`](@ref) for more details. - -# Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the magnitude of - surface tension forces, facilitating the fine-tuning of how surface tension phenomena - are represented in the simulation. -""" -struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension - surface_tension_coefficient::ELTYPE - - function SurfaceTensionAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) - end -end - -@doc raw""" - SurfaceTensionMorris(surface_tension_coefficient=1.0) - -This model implements the surface tension approach described by Morris [Morris2000](@cite). -It calculates surface tension forces based on the curvature of the fluid interface -using particle normals and their divergence, making it suitable for simulating -phenomena like droplet formation and capillary wave dynamics. - -See [`surface_tension`](@ref) for more details. - - -# Keywords -- `surface_tension_coefficient=1.0`: Adjusts the magnitude of the surface tension - forces, enabling tuning of fluid surface behaviors in simulations. -""" -struct SurfaceTensionMorris{ELTYPE} <: AbstractSurfaceTension - surface_tension_coefficient::ELTYPE - - function SurfaceTensionMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) - end -end - -function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles) - return (;) -end - -function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) - curvature = Array{ELTYPE, 1}(undef, nparticles) - return (; curvature) -end - -@doc raw""" - SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0) - -This model implements the momentum-conserving surface tension approach outlined by Morris -[Morris2000](@cite). It calculates surface tension forces using the divergence of a stress -tensor, ensuring exact conservation of linear momentum. This method is particularly -useful for simulations where momentum conservation is critical, though it may require -numerical adjustments at higher resolutions. - -See [`surface_tension`](@ref) for more details. - -# Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the strength of surface tension - forces, allowing fine-tuning to replicate physical behavior. -""" -struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension - surface_tension_coefficient::ELTYPE - - function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) - end -end - -function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, - nparticles) - delta_s = Array{ELTYPE, 1}(undef, nparticles) - # Allocate stress tensor for each particle: NDIMS x NDIMS x nparticles - stress_tensor = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, nparticles) - return (; stress_tensor, delta_s) -end - -@inline function stress_tensor(particle_system::AbstractFluidSystem, particle) - return extract_smatrix(particle_system.cache.stress_tensor, particle_system, particle) -end - -# Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`. -# Currently, specializations reducing this to simple multiplications exist only up -# to a power of three, see -# https://github.com/JuliaLang/julia/blob/34934736fa4dcb30697ac1b23d11d5ad394d6a4d/base/intfuncs.jl#L327-L339 -# By using the `@fastpow` macro, we are consciously trading off some precision in the result -# for enhanced computational speed. This is especially useful in scenarios where performance -# is a higher priority than exact precision. -@fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, distance) - (; surface_tension_coefficient) = surface_tension - - # Eq. 2 - # We only reach this function when `sqrt(eps()) < distance <= support_radius` - if distance > 0.5 * support_radius - # Attractive force - C = (support_radius - distance)^3 * distance^3 - else - # `distance < 0.5 * support_radius` - # Repulsive force - C = 2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64.0 - end - C *= 32.0 / (pi * support_radius^9) - - # Eq. 1 in acceleration form - cohesion_force = -surface_tension_coefficient * m_b * C * pos_diff / distance - - return cohesion_force -end - -@inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, - distance, adhesion_coefficient) - - # The neighborhood search has an `<=` check, but for `distance == support_radius` - # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`. - # TODO Change this in the neighborhood search? - # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19 - distance >= support_radius && return zero(pos_diff) - - distance <= 0.5 * support_radius && return zero(pos_diff) - - # Eq. 7 - A = 0.007 / support_radius^3.25 * - (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25 - - # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as - # `m_b = rho_0 * volume` (Akinci boundary condition treatment) - adhesion_force = -adhesion_coefficient * m_b * A * pos_diff / distance - - return adhesion_force -end - -# Skip -@inline function surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, particle, neighbor, - pos_diff, distance, rho_a, rho_b, grad_kernel) - return zero(pos_diff) -end - -@inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, - surface_tension_b::CohesionForceAkinci, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) - # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) - - m_b = hydrodynamic_mass(neighbor_system, neighbor) - support_radius = compact_support(smoothing_kernel, - smoothing_length(particle_system, particle)) - - return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) -end - -@inline function surface_tension_force(surface_tension_a::SurfaceTensionAkinci, - surface_tension_b::SurfaceTensionAkinci, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, particle, - neighbor, - pos_diff, distance, rho_a, rho_b, grad_kernel) - (; smoothing_kernel) = particle_system - (; surface_tension_coefficient) = surface_tension_a - - smoothing_length_ = smoothing_length(particle_system, particle) - # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) - - m_b = hydrodynamic_mass(neighbor_system, neighbor) - n_a = surface_normal(particle_system, particle) - n_b = surface_normal(neighbor_system, neighbor) - support_radius = compact_support(smoothing_kernel, smoothing_length_) - - return cohesion_force_akinci(surface_tension_a, support_radius, m_b, - pos_diff, distance) .- - (surface_tension_coefficient * (n_a - n_b) * smoothing_length_) -end - -@inline function surface_tension_force(surface_tension_a::SurfaceTensionMorris, - surface_tension_b::SurfaceTensionMorris, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) - (; surface_tension_coefficient) = surface_tension_a - - # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) - - n_a = surface_normal(particle_system, particle) - curvature_a = curvature(particle_system, particle) - - return -surface_tension_coefficient / rho_a * curvature_a * n_a -end - -function compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - return system -end - -# Section 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" -function compute_stress_tensors!(system::AbstractFluidSystem, - ::SurfaceTensionMomentumMorris, - v, u, v_ode, u_ode, semi, t) - (; cache) = system - (; delta_s, stress_tensor) = cache - - # Reset surface stress_tensor - set_zero!(stress_tensor) - - max_delta_s = maximum(delta_s) - NDIMS = ndims(system) - - @trixi_timeit timer() "compute surface stress tensor" begin - @threaded semi for particle in each_integrated_particle(system) - normal = surface_normal(system, particle) - delta_s_particle = delta_s[particle] - if delta_s_particle > eps() - for i in 1:NDIMS, j in 1:NDIMS - delta_ij = (i == j) ? 1 : 0 - stress_tensor[i, j, - particle] = delta_s_particle * - (delta_ij - normal[i] * normal[j]) - - delta_ij * max_delta_s - end - end - end - end - - return system -end - -function compute_surface_delta_function!(system, surface_tension, semi) - return system -end - -# Eq. 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" -function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris, semi) - (; cache) = system - (; delta_s) = cache - - set_zero!(delta_s) - - @threaded semi for particle in each_integrated_particle(system) - delta_s[particle] = norm(surface_normal(system, particle)) - end - return system -end - -@inline function surface_tension_force(surface_tension_a::SurfaceTensionMomentumMorris, - surface_tension_b::SurfaceTensionMomentumMorris, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) - (; surface_tension_coefficient) = surface_tension_a - - # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) - - S_a = stress_tensor(particle_system, particle) - S_b = stress_tensor(neighbor_system, neighbor) - - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - return surface_tension_coefficient * m_b * (S_a + S_b) / (rho_a * rho_b) * grad_kernel -end - -@inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractBoundarySystem, particle, neighbor, - pos_diff, distance) - (; adhesion_coefficient) = neighbor_system - - # No adhesion with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) - - # No reason to calculate the adhesion force if adhesion coefficient is near zero - abs(adhesion_coefficient) < eps() && return zero(pos_diff) - - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - support_radius = compact_support(particle_system.smoothing_kernel, - smoothing_length(particle_system, particle)) - return adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, - adhesion_coefficient) -end - -@inline function adhesion_force(surface_tension, particle_system, neighbor_system, particle, - neighbor, pos_diff, distance) - return zero(pos_diff) -end diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl deleted file mode 100644 index fcca2383f..000000000 --- a/src/schemes/fluid/viscosity.jl +++ /dev/null @@ -1,449 +0,0 @@ -# Unpack the neighboring systems viscosity to dispatch on the viscosity type. -# This function is only necessary to allow `nothing` as viscosity. -# Otherwise, we could just apply the viscosity as a function directly. -@propagate_inbounds function dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - viscosity = viscosity_model(particle_system, neighbor_system) - - return dv_viscosity(viscosity, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) -end - -@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - return viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) -end - -@inline function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - return zero(pos_diff) -end - -@doc raw""" - ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01) - -Artificial viscosity by Monaghan ([Monaghan1992](@cite), [Monaghan1989](@cite)). - -See [`Viscosity`](@ref viscosity_sph) for an overview and comparison of implemented viscosity models. - -# Keywords -- `alpha`: A value of `0.02` is usually used for most simulations. For a relation with the - kinematic viscosity, see description above. -- `beta=0.0`: A value of `0.0` works well for most fluid simulations and simulations - with shocks of moderate strength. In simulations where the Mach number can be - very high, eg. astrophysical calculation, good results can be obtained by - choosing a value of `beta=2.0` and `alpha=1.0`. -- `epsilon=0.01`: Parameter to prevent singularities. -""" -struct ArtificialViscosityMonaghan{ELTYPE} - alpha :: ELTYPE - beta :: ELTYPE - epsilon :: ELTYPE - - function ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01) - new{typeof(alpha)}(alpha, beta, epsilon) - end -end - -@doc raw""" - ViscosityMorris(; nu, epsilon=0.01) - -Viscosity by [Morris (1997)](@cite Morris1997) also used by [Fourtakas (2019)](@cite Fourtakas2019). - -See [`Viscosity`](@ref viscosity_sph) for an overview and comparison of implemented viscosity models. - -# Keywords -- `nu`: Kinematic viscosity -- `epsilon=0.01`: Parameter to prevent singularities -""" -struct ViscosityMorris{ELTYPE} - nu::ELTYPE - epsilon::ELTYPE - - function ViscosityMorris(; nu, epsilon=0.01) - new{typeof(nu)}(nu, epsilon) - end -end - -function kinematic_viscosity(system, viscosity::ViscosityMorris, smoothing_length, - sound_speed) - return viscosity.nu -end - -@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, - ViscosityMorris})(particle_system, - neighbor_system, - v_particle_system, - v_neighbor_system, - particle, neighbor, - pos_diff, distance, - sound_speed, - m_a, m_b, rho_a, rho_b, - grad_kernel) - rho_mean = (rho_a + rho_b) / 2 - - v_a = viscous_velocity(v_particle_system, particle_system, particle) - v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) - v_diff = v_a - v_b - - smoothing_length_particle = smoothing_length(particle_system, particle) - smoothing_length_neighbor = smoothing_length(particle_system, neighbor) - smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - - nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system), - smoothing_length_particle, sound_speed) - nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system), - smoothing_length_neighbor, sound_speed) - - pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length_average, grad_kernel, nu_a, nu_b) - - return m_b * pi_ab -end - -@inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance, - rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b) - (; alpha, beta, epsilon) = viscosity - - # v_ab ⋅ r_ab - vr = dot(v_diff, pos_diff) - - # Monaghan 2005 p. 1741 (doi: 10.1088/0034-4885/68/8/r01): - # "In the case of shock tube problems, it is usual to turn the viscosity on for - # approaching particles and turn it off for receding particles. In this way, the - # viscosity is used for shocks and not rarefactions." - if vr < 0 - mu = h * vr / (distance^2 + epsilon * h^2) - return (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel - end - - return zero(v_diff) -end - -@inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean, - rho_a, rho_b, h, grad_kernel, nu_a, - nu_b) - epsilon = viscosity.epsilon - - mu_a = nu_a * rho_a - mu_b = nu_b * rho_b - - return (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) / - (distance^2 + epsilon * h^2) * v_diff -end - -# See, e.g., -# Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". -# In: Reports on Progress in Physics (2005), pages 1703-1759. -# [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01) -function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan, - smoothing_length, sound_speed) - (; alpha) = viscosity - - return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4) -end - -@doc raw""" - ViscosityAdami(; nu, epsilon=0.01) - -Viscosity by [Adami (2012)](@cite Adami2012). - -See [`Viscosity`](@ref viscosity_sph) for an overview and comparison of implemented viscosity models. - -# Keywords -- `nu`: Kinematic viscosity -- `epsilon=0.01`: Parameter to prevent singularities -""" -struct ViscosityAdami{ELTYPE} - nu::ELTYPE - epsilon::ELTYPE - - function ViscosityAdami(; nu, epsilon=0.01) - new{typeof(nu)}(nu, epsilon) - end -end - -function adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, - m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) - eta_a = nu_a * rho_a - eta_b = nu_b * rho_b - - eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) - - volume_a = m_a / rho_a - volume_b = m_b / rho_b - - # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 - # They argued that the formulation is more flexible because of the possibility to formulate - # different inter-particle averages or to assume different inter-particle distributions. - # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. - # - # TODO: Is there a better formulation to discretize the Laplace operator? - # Because when using this formulation for the pressure acceleration, it is not - # energy conserving. - # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 - visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a - - return visc .* v_diff -end - -@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) - epsilon = viscosity.epsilon - - smoothing_length_particle = smoothing_length(particle_system, particle) - smoothing_length_neighbor = smoothing_length(particle_system, neighbor) - smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - - nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system), - smoothing_length_particle, sound_speed) - nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system), - smoothing_length_neighbor, sound_speed) - - v_a = viscous_velocity(v_particle_system, particle_system, particle) - v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) - v_diff = v_a - v_b - - return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, - m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) -end - -function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, - sound_speed) - return viscosity.nu -end - -@propagate_inbounds function viscous_velocity(v, system, particle) - return current_velocity(v, system, particle) -end - -@doc raw""" - ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) - -Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) -by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. -The effective kinematic viscosity is defined as - -```math -\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, -``` - -with - -```math -\nu_{\text{SGS}} = (C_S h)^2 |S|, -``` - -and an approximation for the strain rate magnitude given by - -```math -|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, -``` - -where: -- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), -- ``h`` is the local smoothing length. - -The effective dynamic viscosities are then computed as -```math -\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, -``` -and averaged as -```math -\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. -``` - -This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. - -# Keywords -- `nu`: Standard kinematic viscosity. -- `C_S`: Smagorinsky constant. -- `epsilon=0.01`: Parameter to prevent singularities -""" -struct ViscosityAdamiSGS{ELTYPE} - nu :: ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] - C_S :: ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] - epsilon :: ELTYPE # Epsilon for singularity prevention [e.g., 0.001] -end - -ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) - -@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, - neighbor_system, - v_particle_system, - v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) - epsilon = viscosity.epsilon - - smoothing_length_particle = smoothing_length(particle_system, particle) - smoothing_length_neighbor = smoothing_length(particle_system, neighbor) - smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - - nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system), - smoothing_length_particle, sound_speed) - nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system), - smoothing_length_neighbor, sound_speed) - - v_a = viscous_velocity(v_particle_system, particle_system, particle) - v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) - v_diff = v_a - v_b - - # ------------------------------------------------------------------------------ - # SGS part: Compute the subgrid-scale eddy viscosity. - # ------------------------------------------------------------------------------ - # In classical LES [Lilly (1967)](@cite Lilly1967) the Smagorinsky model defines: - # ν_SGS = (C_S Δ)^2 |S|, - # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). - # - # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. - # A common low-order surrogate is to approximate the strain‐rate magnitude by a - # finite difference along each particle pair: - # - # |S| ≈ ‖v_ab‖ / (‖r_ab‖ + δ), - # - # where δ regularizes the denominator to avoid singularities when particles are very close. - # - # This yields: - # S_mag = norm(v_diff) / (distance + ε) - # - # and then the Smagorinsky eddy viscosity: - # ν_SGS = (C_S * h̄)^2 * S_mag. - # - S_mag = norm(v_diff) / (distance + epsilon) - nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag - - # Effective kinematic viscosity is the sum of the standard and SGS parts. - nu_a = nu_a + nu_SGS - nu_b = nu_b + nu_SGS - - return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, - m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) -end - -function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, - sound_speed) - return viscosity.nu -end - -@doc raw""" - ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) - -Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997), -by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. -The effective kinematic viscosity is defined as - -```math -\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, -``` - -with - -```math -\nu_{\text{SGS}} = (C_S h)^2 |S|, -``` - -and an approximation for the strain rate magnitude given by - -```math -|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, -``` - -where: -- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), -- ``h`` is the local smoothing length. - -The effective dynamic viscosities are then computed as -```math -\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, -``` -and averaged as -```math -\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. -``` - -This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. - -# Keywords -- `nu`: Standard kinematic viscosity. -- `C_S`: Smagorinsky constant. -- `epsilon=0.01`: Parameter to prevent singularities -""" -struct ViscosityMorrisSGS{ELTYPE} - nu::ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] - C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] - epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] -end - -ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon) - -@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system, - neighbor_system, - v_particle_system, - v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, - m_b, rho_a, rho_b, grad_kernel) - epsilon = viscosity.epsilon - - smoothing_length_particle = smoothing_length(particle_system, particle) - smoothing_length_neighbor = smoothing_length(particle_system, neighbor) - smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - - nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system), - smoothing_length_particle, sound_speed) - nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system), - smoothing_length_neighbor, sound_speed) - - v_a = viscous_velocity(v_particle_system, particle_system, particle) - v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) - v_diff = v_a - v_b - - # SGS part: Compute the subgrid-scale eddy viscosity. - # See comments above for `ViscosityAdamiSGS`. - S_mag = norm(v_diff) / (distance + epsilon) - nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag - - # Effective viscosities include the SGS term. - nu_a_eff = nu_a + nu_SGS - nu_b_eff = nu_b + nu_SGS - - # For the Morris model, dynamic viscosities are: - mu_a = nu_a_eff * rho_a - mu_b = nu_b_eff * rho_b - - force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / - (distance^2 + epsilon * smoothing_length_average^2) * v_diff - return m_b * force_Morris -end - -function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length, - sound_speed) - return viscosity.nu -end diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl deleted file mode 100644 index d1852f78d..000000000 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ /dev/null @@ -1,244 +0,0 @@ -@doc raw""" - AbstractDensityDiffusion - -An abstract supertype of all density diffusion formulations. - -Currently, the following formulations are available: - -| Formulation | Suitable for Steady-State Simulations | Low Computational Cost | -| :------------------------------------------ | :------------------------------------ | :--------------------- | -| [`DensityDiffusionMolteniColagrossi`](@ref) | ❌ | ✅ | -| [`DensityDiffusionFerrari`](@ref) | ❌ | ✅ | -| [`DensityDiffusionAntuono`](@ref) | ✅ | ❌ | - -See [Density Diffusion](@ref density_diffusion) for a comparison and more details. -""" -abstract type AbstractDensityDiffusion end - -@inline density_diffusion(system) = nothing - -# Most density diffusion formulations don't need updating -function update!(density_diffusion, v, u, system, semi) - return density_diffusion -end - -@doc raw""" - DensityDiffusionMolteniColagrossi(; delta) - -The commonly used density diffusion term by [Molteni (2009)](@cite Molteni2009). - -The term ``\psi_{ab}`` in the continuity equation in -[`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) is defined by -```math -\psi_{ab} = 2(\rho_a - \rho_b) \frac{r_{ab}}{\Vert r_{ab} \Vert^2}, -``` -where ``\rho_a`` and ``\rho_b`` denote the densities of particles ``a`` and ``b`` respectively -and ``r_{ab} = r_a - r_b`` is the difference of the coordinates of particles ``a`` and ``b``. - -See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) -for an overview and comparison of implemented density diffusion terms. -""" -struct DensityDiffusionMolteniColagrossi{ELTYPE} <: AbstractDensityDiffusion - delta::ELTYPE - - function DensityDiffusionMolteniColagrossi(; delta) - new{typeof(delta)}(delta) - end -end - -@inline function density_diffusion_psi(::DensityDiffusionMolteniColagrossi, rho_a, rho_b, - pos_diff, distance, system, particle, neighbor) - return 2 * (rho_a - rho_b) * pos_diff / distance^2 -end - -@doc raw""" - DensityDiffusionFerrari() - -A density diffusion term by [Ferrari (2009)](@cite Ferrari2009). - -The term ``\psi_{ab}`` in the continuity equation in -[`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) is defined by -```math -\psi_{ab} = \frac{\rho_a - \rho_b}{h_a + h_b} \frac{r_{ab}}{\Vert r_{ab} \Vert}, -``` -where ``\rho_a`` and ``\rho_b`` denote the densities of particles ``a`` and ``b`` respectively, -``r_{ab} = r_a - r_b`` is the difference of the coordinates of particles ``a`` and ``b`` and -``h_a`` and ``h_b`` are the smoothing lengths of particles ``a`` and ``b`` respectively. - -See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) -for an overview and comparison of implemented density diffusion terms. -""" -struct DensityDiffusionFerrari <: AbstractDensityDiffusion - delta::Int - - # δ is always 1 in this formulation - DensityDiffusionFerrari() = new(1) -end - -@inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b, - pos_diff, distance, system, particle, neighbor) - return ((rho_a - rho_b) / - (smoothing_length(system, particle) + smoothing_length(system, neighbor))) * - pos_diff / distance -end - -@doc raw""" - DensityDiffusionAntuono(initial_condition; delta) - -The commonly used density diffusion terms by [Antuono (2010)](@cite Antuono2010), also referred to as -δ-SPH. The density diffusion term by [Molteni (2009)](@cite Molteni2009) is extended by a second -term, which is nicely written down by [Antuono (2012)](@cite Antuono2012). - -The term ``\psi_{ab}`` in the continuity equation in -[`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) is defined by -```math -\psi_{ab} = 2\left(\rho_a - \rho_b - \frac{1}{2}\big(\nabla\rho^L_a + \nabla\rho^L_b\big) \cdot r_{ab}\right) - \frac{r_{ab}}{\Vert r_{ab} \Vert^2}, -``` -where ``\rho_a`` and ``\rho_b`` denote the densities of particles ``a`` and ``b`` respectively -and ``r_{ab} = r_a - r_b`` is the difference of the coordinates of particles ``a`` and ``b``. -The symbol ``\nabla\rho^L_a`` denotes the renormalized density gradient defined as -```math -\nabla\rho^L_a = -\sum_b (\rho_a - \rho_b) V_b L_a \nabla W_{ab} -``` -with -```math -L_a := \left( -\sum_{b} V_b r_{ab} \otimes \nabla W_{ab} \right)^{-1} \in \R^{d \times d}, -``` -where ``d`` is the number of dimensions. - -See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) -for an overview and comparison of implemented density diffusion terms. -""" -struct DensityDiffusionAntuono{NDIMS, ELTYPE, ARRAY2D, ARRAY3D} <: AbstractDensityDiffusion - delta :: ELTYPE - correction_matrix :: ARRAY3D # Array{ELTYPE, 3}: [i, j, particle] - normalized_density_gradient :: ARRAY2D # Array{ELTYPE, 2}: [i, particle] - - function DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient) - new{size(correction_matrix, 1), typeof(delta), - typeof(normalized_density_gradient), - typeof(correction_matrix)}(delta, correction_matrix, - normalized_density_gradient) - end -end - -function DensityDiffusionAntuono(initial_condition; delta) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, - nparticles(initial_condition)) - - normalized_density_gradient = Array{ELTYPE, 2}(undef, NDIMS, - nparticles(initial_condition)) - - return DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient) -end - -@inline Base.ndims(::DensityDiffusionAntuono{NDIMS}) where {NDIMS} = NDIMS - -function Base.show(io::IO, density_diffusion::DensityDiffusionAntuono) - @nospecialize density_diffusion # reduce precompilation time - - print(io, "DensityDiffusionAntuono(") - print(io, density_diffusion.delta) - print(io, ")") -end - -function allocate_buffer(initial_condition, density_diffusion, buffer) - return allocate_buffer(initial_condition, buffer), density_diffusion -end - -function allocate_buffer(ic, dd::DensityDiffusionAntuono, buffer::SystemBuffer) - initial_condition = allocate_buffer(ic, buffer) - return initial_condition, DensityDiffusionAntuono(initial_condition; delta=dd.delta) -end - -@inline function density_diffusion_psi(density_diffusion::DensityDiffusionAntuono, - rho_a, rho_b, pos_diff, distance, system, - particle, neighbor) - (; normalized_density_gradient) = density_diffusion - - normalized_gradient_a = extract_svector(normalized_density_gradient, system, particle) - normalized_gradient_b = extract_svector(normalized_density_gradient, system, neighbor) - - # Fist term by Molteni & Colagrossi - result = 2 * (rho_a - rho_b) - - # Second correction term - result -= dot(normalized_gradient_a + normalized_gradient_b, pos_diff) - - return result * pos_diff / distance^2 -end - -function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi) - (; normalized_density_gradient) = density_diffusion - - # Compute correction matrix - density_fun = @inline(particle->current_density(v, system, particle)) - system_coords = current_coordinates(u, system) - - compute_gradient_correction_matrix!(density_diffusion.correction_matrix, - system, system_coords, density_fun, semi) - - # Compute normalized density gradient - set_zero!(normalized_density_gradient) - - foreach_point_neighbor(system, system, system_coords, system_coords, semi; - points=each_integrated_particle(system)) do particle, neighbor, - pos_diff, distance - # Density diffusion terms are all zero for distance zero. - # See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(system)^2) && return - - rho_a = current_density(v, system, particle) - rho_b = current_density(v, system, neighbor) - - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - L = correction_matrix(density_diffusion, particle) - - m_b = hydrodynamic_mass(system, neighbor) - volume_b = m_b / rho_b - - normalized_gradient = -(rho_a - rho_b) * L * grad_kernel * volume_b - - for i in eachindex(normalized_gradient) - normalized_density_gradient[i, particle] += normalized_gradient[i] - end - end - - return density_diffusion -end - -@propagate_inbounds function density_diffusion!(dv, - density_diffusion::AbstractDensityDiffusion, - v_particle_system, particle, neighbor, - pos_diff, distance, m_b, rho_a, rho_b, - particle_system::Union{AbstractFluidSystem, - OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - grad_kernel) - # Density diffusion terms are all zero for distance zero. - # See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return - - (; delta) = density_diffusion - sound_speed = system_sound_speed(particle_system) - - volume_b = m_b / rho_b - - psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance, - particle_system, particle, neighbor) - density_diffusion_term = dot(psi, grad_kernel) * volume_b - - smoothing_length_avg = (smoothing_length(particle_system, particle) + - smoothing_length(particle_system, neighbor)) / 2 - dv[end, particle] += delta * smoothing_length_avg * sound_speed * - density_diffusion_term -end - -# Density diffusion `nothing` or interaction other than fluid-fluid -@inline function density_diffusion!(dv, density_diffusion, v_particle_system, particle, - neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, grad_kernel) - return dv -end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl deleted file mode 100644 index 5462a663a..000000000 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ /dev/null @@ -1,442 +0,0 @@ -""" - WeaklyCompressibleSPHSystem(initial_condition, - density_calculator, state_equation, - smoothing_kernel, smoothing_length; - acceleration=ntuple(_ -> 0.0, NDIMS), - viscosity=nothing, density_diffusion=nothing, - pressure_acceleration=nothing, - shifting_technique=nothing, - buffer_size=nothing, - correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, - reference_particle_spacing=0.0)) - -System for particles of a fluid. -The weakly compressible SPH (WCSPH) scheme is used, wherein a stiff equation of state -generates large pressure changes for small density variations. -See [Weakly Compressible SPH](@ref wcsph) for more details on the method. - -# Arguments -- `initial_condition`: [`InitialCondition`](@ref) representing the system's particles. -- `density_calculator`: Density calculator for the system. - See [`ContinuityDensity`](@ref) and [`SummationDensity`](@ref). -- `state_equation`: Equation of state for the system. See [`StateEquationCole`](@ref). -- `smoothing_kernel`: Smoothing kernel to be used for this system. - See [Smoothing Kernels](@ref smoothing_kernel). -- `smoothing_length`: Smoothing length to be used for this system. - See [Smoothing Kernels](@ref smoothing_kernel). - -# Keywords -- `acceleration`: Acceleration vector for the system. (default: zero vector) -- `viscosity`: Viscosity model for this system (default: no viscosity). - See [`ArtificialViscosityMonaghan`](@ref) or [`ViscosityAdami`](@ref). -- `density_diffusion`: Density diffusion terms for this system. - See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion). -- `pressure_acceleration`: Pressure acceleration formulation for this system. - By default, the correct formulation is chosen based on the - density calculator and the correction method. - To use [Tensile Instability Control](@ref tic), pass - [`tensile_instability_control`](@ref) here. -- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity - formulation](@ref transport_velocity_formulation) to use - with this system. Default is no shifting. -- `buffer_size`: Number of buffer particles. - This is needed when simulating with [`OpenBoundarySystem`](@ref). -- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) -- `source_terms`: Additional source terms for this system. Has to be either `nothing` - (by default), or a function of `(coords, velocity, density, pressure, t)` - (which are the quantities of a single particle), returning a `Tuple` - or `SVector` that is to be added to the acceleration of that particle. - See, for example, [`SourceTermDamping`](@ref). - Note that these source terms will not be used in the calculation of the - boundary pressure when using a boundary with - [`BoundaryModelDummyParticles`](@ref) and [`AdamiPressureExtrapolation`](@ref). - The keyword argument `acceleration` should be used instead for - gravity-like source terms. -- `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) -- `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: The value used to initialize the color of particles in the system. -""" -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, - PF, SC, ST, B, SRFT, SRFN, PR, - C} <: AbstractFluidSystem{NDIMS} - initial_condition :: IC - mass :: MA # Array{ELTYPE, 1} - pressure :: P # Array{ELTYPE, 1} - density_calculator :: DC - state_equation :: SE - smoothing_kernel :: K - acceleration :: SVector{NDIMS, ELTYPE} - viscosity :: V - density_diffusion :: DD - correction :: COR - pressure_acceleration_formulation :: PF - shifting_technique :: SC - source_terms :: ST - surface_tension :: SRFT - surface_normal_method :: SRFN - buffer :: B - particle_refinement :: PR # TODO - cache :: C -end - -# The default constructor needs to be accessible for Adapt.jl to work with this struct. -# See the comments in general/gpu.jl for more details. -function WeaklyCompressibleSPHSystem(initial_condition, - density_calculator, state_equation, - smoothing_kernel, smoothing_length; - acceleration=ntuple(_ -> zero(eltype(initial_condition)), - ndims(smoothing_kernel)), - viscosity=nothing, density_diffusion=nothing, - pressure_acceleration=nothing, - shifting_technique=nothing, - buffer_size=nothing, - correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, - reference_particle_spacing=0, color_value=1) - buffer = isnothing(buffer_size) ? nothing : - SystemBuffer(nparticles(initial_condition), buffer_size) - - particle_refinement = nothing # TODO - - initial_condition, - density_diffusion = allocate_buffer(initial_condition, - density_diffusion, buffer) - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - n_particles = nparticles(initial_condition) - - mass = copy(initial_condition.mass) - pressure = similar(initial_condition.pressure) - - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end - - # Make acceleration an SVector - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end - - if correction isa ShepardKernelCorrection && - density_calculator isa ContinuityDensity - throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) - end - - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end - - if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) - end - - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - correction) - - cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_correction(correction, initial_condition.density, NDIMS, - n_particles)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., - create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, - n_particles)..., - create_cache_refinement(initial_condition, particle_refinement, - smoothing_length)..., - create_cache_shifting(initial_condition, shifting_technique)..., - color=Int(color_value)) - - # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` - if reference_particle_spacing > 0 - # `reference_particle_spacing` has to be set for surface normals to be determined - cache = (; - cache..., # Existing cache fields - reference_particle_spacing=reference_particle_spacing) - end - - return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, - density_calculator, state_equation, - smoothing_kernel, acceleration_, viscosity, - density_diffusion, correction, pressure_acceleration, - shifting_technique, source_terms, surface_tension, - surface_normal_method, buffer, particle_refinement, - cache) -end - -function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) - @nospecialize system # reduce precompilation time - - print(io, "WeaklyCompressibleSPHSystem{", ndims(system), "}(") - print(io, system.density_calculator) - print(io, ", ", system.correction) - print(io, ", ", system.state_equation) - print(io, ", ", system.smoothing_kernel) - print(io, ", ", system.viscosity) - print(io, ", ", system.density_diffusion) - print(io, ", ", system.shifting_technique) - print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal - print(io, ", ", system.color) - end - print(io, ", ", system.acceleration) - print(io, ", ", system.source_terms) - print(io, ") with ", nparticles(system), " particles") -end - -function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSystem) - @nospecialize system # reduce precompilation time - - if get(io, :compact, false) - show(io, system) - else - summary_header(io, "WeaklyCompressibleSPHSystem{$(ndims(system))}") - if system.buffer isa SystemBuffer - summary_line(io, "#particles", nparticles(system)) - summary_line(io, "#buffer_particles", system.buffer.buffer_size) - else - summary_line(io, "#particles", nparticles(system)) - end - summary_line(io, "density calculator", - system.density_calculator |> typeof |> nameof) - summary_line(io, "correction method", - system.correction |> typeof |> nameof) - summary_line(io, "state equation", system.state_equation |> typeof |> nameof) - summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "viscosity", system.viscosity) - summary_line(io, "density diffusion", system.density_diffusion) - summary_line(io, "shifting technique", system.shifting_technique) - summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal - summary_line(io, "color", system.cache.color) - end - summary_line(io, "acceleration", system.acceleration) - summary_line(io, "source terms", system.source_terms |> typeof |> nameof) - summary_footer(io) - end -end - -@inline function Base.eltype(::WeaklyCompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE} - return ELTYPE -end - -@inline function v_nvariables(system::WeaklyCompressibleSPHSystem) - return v_nvariables(system, system.density_calculator) -end - -@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::SummationDensity) - return ndims(system) -end - -@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity) - return ndims(system) + 1 -end - -@inline buffer(system::WeaklyCompressibleSPHSystem) = system.buffer - -system_correction(system::WeaklyCompressibleSPHSystem) = system.correction - -@inline function current_velocity(v, system::WeaklyCompressibleSPHSystem) - return current_velocity(v, system.density_calculator, system) -end - -@inline function current_velocity(v, ::SummationDensity, - system::WeaklyCompressibleSPHSystem) - # When using `SummationDensity`, `v` contains only the velocity - return v -end - -@inline function current_velocity(v, ::ContinuityDensity, - system::WeaklyCompressibleSPHSystem) - # When using `ContinuityDensity`, the velocity is stored - # in the first `ndims(system)` rows of `v`. - return view(v, 1:ndims(system), :) -end - -@inline function current_density(v, system::WeaklyCompressibleSPHSystem) - return current_density(v, system.density_calculator, system) -end - -@inline function current_density(v, ::SummationDensity, - system::WeaklyCompressibleSPHSystem) - # When using `SummationDensity`, the density is stored in the cache - return system.cache.density -end - -@inline function current_density(v, ::ContinuityDensity, - system::WeaklyCompressibleSPHSystem) - # When using `ContinuityDensity`, the density is stored in the last row of `v` - return view(v, size(v, 1), :) -end - -@inline function current_pressure(v, system::WeaklyCompressibleSPHSystem) - return system.pressure -end - -@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed - -@inline shifting_technique(system::WeaklyCompressibleSPHSystem) = system.shifting_technique - -@inline density_diffusion(system::WeaklyCompressibleSPHSystem) = system.density_diffusion - -function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, - v_ode, u_ode, semi, t) - (; density_calculator, density_diffusion, correction) = system - - compute_density!(system, u, u_ode, semi, density_calculator) - - @trixi_timeit timer() "update density diffusion" update!(density_diffusion, v, u, - system, semi) - - return system -end - -function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_normal_method, surface_tension) = system - - compute_pressure!(system, v, semi) - - # These are only computed when using corrections - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) - compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) - # `kernel_correct_density!` only performed for `SummationDensity` - kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, - density_calculator) - - # These are only computed when using surface tension - compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) - compute_surface_delta_function!(system, surface_tension, semi) - return system -end - -function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; surface_tension) = system - - # Surface normal of neighbor and boundary needs to have been calculated already - compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) -end - -function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, - semi, correction, density_calculator) - return system -end - -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) = 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_kernel) -end - -function reinit_density!(vu_ode, semi) - v_ode, u_ode = vu_ode.x - - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - reinit_density!(system, v, u, v_ode, u_ode, semi) - end - - return vu_ode -end - -function reinit_density!(system::WeaklyCompressibleSPHSystem, v, u, - v_ode, u_ode, semi) - # Compute density with `SummationDensity` and store the result in `v`, - # overwriting the previous integrated density. - summation_density!(system, semi, u, u_ode, v[end, :]) - - # Apply `ShepardKernelCorrection` - kernel_correction_coefficient = zeros(size(v[end, :])) - compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi, - kernel_correction_coefficient) - v[end, :] ./= kernel_correction_coefficient - - compute_pressure!(system, v, semi) - - return system -end - -function reinit_density!(system, v, u, v_ode, u_ode, semi) - return system -end - -function compute_pressure!(system, v, semi) - @threaded semi for particle in eachparticle(system) - apply_state_equation!(system, current_density(v, system, particle), 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. -@inline function apply_state_equation!(system::WeaklyCompressibleSPHSystem, density, - particle) - system.pressure[particle] = system.state_equation(density) -end - -function write_v0!(v0, system::WeaklyCompressibleSPHSystem, ::ContinuityDensity) - # Note that `.=` is very slightly faster, but not GPU-compatible - v0[end, :] = system.initial_condition.density - - return v0 -end - -function restart_with!(system::WeaklyCompressibleSPHSystem, v, u) - for particle in each_integrated_particle(system) - system.initial_condition.coordinates[:, particle] .= u[:, particle] - system.initial_condition.velocity[:, particle] .= v[1:ndims(system), particle] - end - - restart_with!(system, system.density_calculator, v, u) -end - -function restart_with!(system, ::SummationDensity, v, u) - return system -end - -function restart_with!(system, ::ContinuityDensity, v, u) - for particle in each_integrated_particle(system) - system.initial_condition.density[particle] = v[end, particle] - end - - return system -end - -@inline function correction_matrix(system::WeaklyCompressibleSPHSystem, particle) - extract_smatrix(system.cache.correction_matrix, system, particle) -end - -@inline function curvature(particle_system::AbstractFluidSystem, particle) - (; cache) = particle_system - return cache.curvature[particle] -end